{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Modeling and Simulation in Python\n",
    "\n",
    "Case study.\n",
    "\n",
    "Copyright 2017 Allen Downey\n",
    "\n",
    "License: [Creative Commons Attribution 4.0 International](https://creativecommons.org/licenses/by/4.0)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Configure Jupyter so figures appear in the notebook\n",
    "%matplotlib inline\n",
    "\n",
    "# Configure Jupyter to display the assigned value after an assignment\n",
    "%config InteractiveShell.ast_node_interactivity='last_expr_or_assign'\n",
    "\n",
    "# import functions from the modsim.py module\n",
    "from modsim import *"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Electric car"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Olin Electric Motorsports](https://www.olinelectricmotorsports.com/) is a club at Olin College that designs and builds electric cars, and participates in the [Formula SAE Electric](https://www.sae.org/attend/student-events/formula-sae-electric) competition.\n",
    "\n",
    "The goal of this case study is to use simulation to guide the design of a car intended to accelerate from standing to 100 kph as quickly as possible.  The [world record for this event](https://www.youtube.com/watch?annotation_id=annotation_2297602723&feature=iv&src_vid=I-NCH8ct24U&v=n2XiCYA3C9s), using a car that meets the competition requirements, is 1.513 seconds.\n",
    "\n",
    "We'll start with a simple model that takes into account the characteristics of the motor and vehicle:\n",
    "\n",
    "* The motor is an [Emrax 228 high voltage axial flux synchronous permanent magnet motor](http://emrax.com/products/emrax-228/); according to the [data sheet](http://emrax.com/wp-content/uploads/2017/01/emrax_228_technical_data_4.5.pdf), its maximum torque is 240 Nm, at 0 rpm.  But maximum torque decreases with motor speed; at 5000 rpm, maximum torque is 216 Nm.\n",
    "\n",
    "* The motor is connected to the drive axle with a chain drive with speed ratio 13:60 or 1:4.6; that is, the axle rotates once for each 4.6 rotations of the motor.\n",
    "\n",
    "* The radius of the tires is 0.26 meters.\n",
    "\n",
    "* The weight of the vehicle, including driver, is 300 kg.\n",
    "\n",
    "To start, we will assume no slipping between the tires and the road surface, no air resistance, and no rolling resistance.  Then we will relax these assumptions one at a time.\n",
    "\n",
    "* First we'll add drag, assuming that the frontal area of the vehicle is 0.6 square meters, with coefficient of drag 0.6.\n",
    "\n",
    "* Next we'll add rolling resistance, assuming a coefficient of 0.2.\n",
    "\n",
    "* Finally we'll compute the peak acceleration to see if the \"no slip\" assumption is credible.\n",
    "\n",
    "We'll use this model to estimate the potential benefit of possible design improvements, including decreasing drag and rolling resistance, or increasing the speed ratio.\n",
    "\n",
    "I'll start by loading the units we need."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "revolutions_per_minute"
      ],
      "text/latex": [
       "$\\mathrm{revolutions_per_minute}$"
      ],
      "text/plain": [
       "<Unit('revolutions_per_minute')>"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "radian = UNITS.radian\n",
    "m = UNITS.meter\n",
    "s = UNITS.second\n",
    "minute = UNITS.minute\n",
    "hour = UNITS.hour\n",
    "km = UNITS.kilometer\n",
    "kg = UNITS.kilogram\n",
    "N = UNITS.newton\n",
    "rpm = UNITS.rpm"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And store the parameters in a `Params` object."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>values</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>r_wheel</th>\n",
       "      <td>0.26 meter</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>speed_ratio</th>\n",
       "      <td>0.216667</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>C_rr</th>\n",
       "      <td>0.2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>C_d</th>\n",
       "      <td>0.5</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>area</th>\n",
       "      <td>0.6 meter ** 2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>rho</th>\n",
       "      <td>1.2 kilogram / meter ** 3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>mass</th>\n",
       "      <td>300 kilogram</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "r_wheel                       0.26 meter\n",
       "speed_ratio                     0.216667\n",
       "C_rr                                 0.2\n",
       "C_d                                  0.5\n",
       "area                      0.6 meter ** 2\n",
       "rho            1.2 kilogram / meter ** 3\n",
       "mass                        300 kilogram\n",
       "dtype: object"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "params = Params(r_wheel=0.26 * m,\n",
    "                speed_ratio=13/60,\n",
    "                C_rr=0.2,\n",
    "                C_d=0.5,\n",
    "                area=0.6*m**2,\n",
    "                rho=1.2*kg/m**3,\n",
    "                mass=300*kg)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`make_system` creates the initial state, `init`, and constructs an `interp1d` object that represents torque as a function of motor speed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def make_system(params):\n",
    "    \"\"\"Make a system object.\n",
    "    \n",
    "    params: Params object\n",
    "    \n",
    "    returns: System object\n",
    "    \"\"\"\n",
    "    init = State(x=0*m, v=0*m/s)\n",
    "    \n",
    "    rpms = [0, 2000, 5000]\n",
    "    torques = [240, 240, 216]\n",
    "    interpolate_torque = interpolate(Series(torques, rpms))\n",
    "    \n",
    "    return System(params, init=init,\n",
    "                  interpolate_torque=interpolate_torque,\n",
    "                  t_end=3*s)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Testing `make_system`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>values</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>r_wheel</th>\n",
       "      <td>0.26 meter</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>speed_ratio</th>\n",
       "      <td>0.216667</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>C_rr</th>\n",
       "      <td>0.2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>C_d</th>\n",
       "      <td>0.5</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>area</th>\n",
       "      <td>0.6 meter ** 2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>rho</th>\n",
       "      <td>1.2 kilogram / meter ** 3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>mass</th>\n",
       "      <td>300 kilogram</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>init</th>\n",
       "      <td>x               0 meter\n",
       "v    0.0 meter / secon...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>interpolate_torque</th>\n",
       "      <td>&lt;function interpolate.&lt;locals&gt;.wrapper at 0x7f...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>unit_torque</th>\n",
       "      <td>1.0 meter * newton</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>t_end</th>\n",
       "      <td>3 second</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "r_wheel                                                      0.26 meter\n",
       "speed_ratio                                                    0.216667\n",
       "C_rr                                                                0.2\n",
       "C_d                                                                 0.5\n",
       "area                                                     0.6 meter ** 2\n",
       "rho                                           1.2 kilogram / meter ** 3\n",
       "mass                                                       300 kilogram\n",
       "init                  x               0 meter\n",
       "v    0.0 meter / secon...\n",
       "interpolate_torque    <function interpolate.<locals>.wrapper at 0x7f...\n",
       "unit_torque                                          1.0 meter * newton\n",
       "t_end                                                          3 second\n",
       "dtype: object"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "system = make_system(params)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>values</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>x</th>\n",
       "      <td>0 meter</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>v</th>\n",
       "      <td>0.0 meter / second</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "x               0 meter\n",
       "v    0.0 meter / second\n",
       "dtype: object"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "system.init"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Torque and speed\n",
    "\n",
    "The relationship between torque and motor speed is taken from the [Emrax 228 data sheet](http://emrax.com/wp-content/uploads/2017/01/emrax_228_technical_data_4.5.pdf).  The following functions reproduce the red dotted line that represents peak torque, which can only be sustained for a few seconds before the motor overheats."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "def compute_torque(omega, system):\n",
    "    \"\"\"Maximum peak torque as a function of motor speed.\n",
    "    \n",
    "    omega: motor speed in radian/s\n",
    "    system: System object\n",
    "    \n",
    "    returns: torque in Nm\n",
    "    \"\"\"\n",
    "    factor = (1 * radian / s).to(rpm)\n",
    "    x = magnitude(omega * factor)\n",
    "    return system.interpolate_torque(x) * N * m"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "240.0 meter newton"
      ],
      "text/latex": [
       "$240.0\\ \\mathrm{meter} \\cdot \\mathrm{newton}$"
      ],
      "text/plain": [
       "240.0 <Unit('meter * newton')>"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "compute_torque(0*radian/s, system)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "216.0 meter newton"
      ],
      "text/latex": [
       "$216.0\\ \\mathrm{meter} \\cdot \\mathrm{newton}$"
      ],
      "text/plain": [
       "216.0 <Unit('meter * newton')>"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "omega = (5000 * rpm).to(radian/s)\n",
    "compute_torque(omega, system)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot the whole curve."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO3dd5hU5fnG8e/sbKNX6U0QHgFBelnsRo2x/IwtdkQFVklTEyuCFMXEmJgYDWBHo0RjwRKDxkQjHRSwIA+ogCAiIL3sLlt+f5whGTewDGX27M7en+vai5lT70PIPr7nvOd9IyUlJYiIiFQ0aWEHEBER2RMVKBERqZBUoEREpEJSgRIRkQopPewAyWZmWUBv4GugKOQ4IiLyXVGgKTDX3fPjV6R8gSIoTu+FHUJERMp0LDAtfkFVKFBfA/z5z3+mSZMmYWcREZE4a9as4dJLL4XY7+p4VaFAFQE0adKEFi1ahJ1FRET27H8ewaiThIiIVEgqUCIiUiGpQImISIWkAiUiIhVSuXaSMLNTgHuA9sBa4F53n1BqmweALu5+QtyyC4G7CfrKvwtc6e5ryyu3iIiUv3JrQZlZS+AFYCxQF7gYGGdmp8Vtczpwban9OgGPAlcCDYClwOTySS0iImEpzxZUG+AZd38p9n2umb0DDACmmtlhwAPAn4AucftdBrzq7tMAzOxWYKOZtXf3pckOvTO/kE+XbaBY05KkhLRIhPT0COnRtLifuO/p/7ssLS0SdmyRKqncCpS7v0fciA5mVp/gzeGnYoseI7j9V5PvFqhOwLy44+wws5WxbZJeoB54bgHvLfgq2aeRCiwtLUJ6WoT09DTq1szi2vO60q1Do7BjiaS8UF7UNbM6wCvAbGCKmV0HFLv7I2b281Kb1wR2lFq2A6ie/KRwYs8W7MwvVAsqRRQXl1BUVEJhUTG7ioopKiqmsKiYwsISCouLKSyMfY9tU1hUTHFxCQXFJRQUFrMjr5Axj81h9JD+dG7bIOzLEUlp5V6gzKwDMAVYBFwKGPBLoM9edtkOVCu1rDqwLVkZ4/Xu1ITenTREUlVVUlJCcXEJhcUl7Cos5rFXPuatOV8y6pFZjM3NoUOremFHFElZ5drN3MyOI2g1vQyc7+55wLlAI2CpmW0CxgHHxD5DUMgs7hjVgVax5SJJFYlEiEbTyMqIUrNaBsMu6MZx3ZuzM7+QkRNnsmz15rAjiqSscmtBmVk74DXgdnd/YPdyd78LuCtuu58D58R1M38GmGZmJwAzCQrYfHdfUk7RRf4jmhbh+ot7kF9QxOxP1nDHhBmMu+4YWjauFXY0kZRTni2oYUAtgq7l2+J+flXWTu7+EXAVMB5YD3QGLkh6WpG9SI+mcfMVveje4TA2byvgjgkzWPPt9rBjiaScSEmKP/w3szbAsrffflujmcshlVdQyJ0Pz+KTL76lUf3q/GrYMTSsW/pxqYiUZdWqVZx88skAh7v78vh1GupI5ABlZ6Yz4uq+dGhVl7UbdjB8/HQ2bs0LO5ZIylCBEjkI1bMzuHNwf9o0rc1X67YzYsJMtu4oCDuWSEpQgRI5SLWqZzJmaA4tGtVk+ddbGDlxJjvydoUdS6TSU4ESOQTq1spibG4OjetXZ+nKTYx6ZBZ5+YVhxxKp1FSgRA6RBnWqMTY3h4Z1slm0bAN3PTGHgl3/M4u1iCRIBUrkEGrSoAZjcnOoWzOLBUvW8atJ8ygsKg47lkilpAIlcoi1aFSLMbk51KyWwZxFa/jtMx9QVJzar3OIJIMKlEgStGlam1FD+lMtK533FnzFA8/Np1hFSmS/qECJJEmHVvUYeU0/sjKjvD13JRNf/ohUfzFe5FBSgRJJos5tGzB8UB/So2m8Pn0ZT76+SEVKJEEqUCJJ1q1DI24d2JtoWoQX/vUZk9/SOMciiVCBEikHfTo34cZLe5IWgWemLubldz8LO5JIhacCJVJOju3WnJ9c2B2AR1/5hDdmLAs5kUjFpgIlUo6+16cVued2BeChFz7kn/O+DDmRSMWlAiVSzs4YcDiDzuwEwO8nz2f6wtUhJxKpmFSgREJw7ontufhUo7gE7n16HnMXrQk7kkiFowIlEpKLTzXOOb4dRcUljHtyLguXrgs7kkiFogIlEpJIJMJVZ3Xm9P5t2FVYzNjHZvPpsg1hxxKpMFSgREIUiUTIPbcrJ/VqSV5BEXc+MpPPVm4KO5ZIhaACJRKytLQIP72wGwOObsaOvEJGTJzBiq+3hB1LJHQqUCIVQDSaxo2X9KRXx8Zs3bGL4RNmsHrdtrBjiYRKBUqkgshIT+PWgb3pekRDNm3N5/bxM1i7YUfYsURCowIlUoFkZkQZflVfOrapz/pNOxk+fgbfbt4ZdiyRUKhAiVQw1bLSGXlNP9q1qMPX327njgkz2LwtP+xYIuVOBUqkAqpRLYNRg/vTqkktVn6zjRETZrJtR0HYsUTKlQqUSAVVp2YWY4fm0LRhDb5YvZk7H5nFjrxdYccSKTf7XaDMLJKMICLyv+rVzmZsbg6H1auGr9jI2MfmkL+rKOxYIuUifV8bmFlNYBBwOtALqG9mxcA6YC7wGvCsu29PZlCRqqpRverclTuAWx58j48+X8/dT8xh+KA+ZKRHw44mklR7bUGZWYaZjQBWAVcA7wM/A34A/B9wC/AZMAT40sxGmFlm8iOLVD1NG9ZgzNAcatfI5IPFa7n36fcpKioOO5ZIUpXVgpoN/APo4u4r97LNUwBm1gG4DpgDdDukCUUEgFZNajN6SH9uHz+DmR99zf2T5/Pzi3sQTdNdd0lNZT2DOsvdbyqjOP2Huy9x958DZxy6aCJSWrsWdblzcD+qZUV554NV/OmFhZSUlIQdSyQp9lqg3P2r/T3YgewjIvvnyNb1ueOqfmSmpzF11goemfKxipSkpH12kgAws9bAKKAzkFV6vbt3PcS5RKQMXY5oyG2D+jD2sdm88t4XZGelc/npHcOOJXJIJVSggGeA+sBkIO9AT2ZmpwD3AO2BtcC97j7BzJoAfwJOAEoIegb+xN03x/abBFwIFMYdrqu7f3GgWUQqu55HNuamy3txz6R5PPePJWRnRrng5A5hxxI5ZBItUN2AHHdfeKAnMrOWwAvAQGAK0BOYambLgZ8QFKxmQEZsuzHAT2O79wDOcfe/H+j5RVJR/y7NuP6i7vz22Q+Y9LdPycqMcvax7cKOJXJIJFqgFgJNY38eqDbAM+7+Uuz7XDN7BxgAnAvg7gVm1hCoAawHMLNqwJHAgoM4t0jKOqFnS/J3FfHH5xfy8Msfk52Zzql9W4cdS+SgJVqgrgammNmzwDLgOy9guPukfR3A3d8D3tv93czqA8cCT7l7QWzZs8CPgEXAQ7FNuxHc2nvYzPoBK4ER7v5agtlFUt5p/dqQX1DEw1M+5o/PLyAzI8oJPVqEHUvkoCRaoK4AjgBuBEqP/V8C7LNAxTOzOsArBO9aTYlbNQi4FngMeBE4DqhFUNhGEbTgzgaeM7P+B3PLUSTVnH1cO/IKinjqjU/53bMfkJURpX+XpmHHEjlgiRaoa4Fcd594sCeMvdQ7haCVdKm7/6c15u55QJ6Z3QQsNbP67v4m8GbcIV4ws0EEhUoFSiTOhd/rQF5BIc+/vZRfPzWPO67qS48jG4UdS+SAJDpYbD7wr4M9mZkdR9Bqehk4P1aQMLO5ZnZi3KZZBLf1tpvZWWY2sNShMjmI3oQiqezy0zty1rFtKSwq5q7HZ/PR5+vDjiRyQBItUHcDd5lZgwM9kZm1I+g+PsLdb3X3+DcL3wdGmVn92Dl+A0xy93wgCvzezPqYWdTMLgFygL8caBaRVBaJRBj8f0dxat/WFBQWM+bRWfiKDWHHEtlvid7iuwzoApxnZluA70xK4+6J3EMYRvA8aZyZjYtb/iDBs637gMUELafngVtjx37ZzG4HngWaxLY5092/TDC7SJUTiUS47vyjyS8o4t35qxj58CzuvnYAbZvXCTuaSMISLVB/PNgTufsNwA1lbJIb+9nTvg8SFDIRSVA0LcLPL+5O/q5CZn28hjsmzOCeYcfQsnGtsKOJJCShAuXuTyY7iIgceunRNG66vBdjH5/DB4vXMnz8dO4ZdixNG9YIO5rIPmnKd5EUl5Ee5daBvTmqXQM2bMln+PjprNtY+m0RkYpHBUqkCsjOTOeOq/pireqxduNOho+fzsYt6ggrFZsKlEgVUT07gzsH96NtszqsXr+dOybMYMv2grBjieyVCpRIFVKzeiajh/anZeOarFizlZETZ7B956597ygSgr12kjCzEQkeo8TdxxyiPCKSZHVqZjFmaA63PDiNz1ZtZtQjsxg9pD/ZWYl26hUpH2X9izxrH/s2I3gvqZhgagwRqSQa1KnG2NwB3PLgND5dvoGxj89mxNX9yMyIhh1N5D/2WqDcvfeelptZOnATMBz4CBiSnGgikkyN61dnbG7Qklq4dD33TJrLrQP7kJGuO/9SMezXv0Qz6w98ANwGjAR6uvucZAQTkeRrflhNxg7NoVb1DOYu+ob7nnmfoqLife8oUg4SKlBmVtvMHgKmAauBLu5+r7sXJTWdiCRd66a1GT0kh+rZ6UxfuJo/PLeA4uKSfe8okmT7LFBmdj7B+HfnApe5+/fdfVnSk4lIuTmiZV1GXtOPrMwo/5y3kvEvfUhJiYqUhGuvBcrMWprZq8BkglHIzd2fLbdkIlKuOh3egDsG9SUjPY03Zizn8dcWqUhJqMrqxbcIqA58CdQlmHJ9jxu6+4WHPpqIlLejOxzGrQN7c9fjc3jpnc/IzoxyyWlHhh1LqqiybvG9QDCV+zvA9n38iEiK6N2pCb+4rCdpEXj2TefFfy0NO5JUUWV1M7+yHHOISAVyzNHNKdhVxO+enc/jry0iKzOdMwYcHnYsqWLKegY10syqJXogM6tlZqMPTSwRCdtJvVpx3XldARj/4of8Y47mCJXyVdYzqC3AJ2b2V+BFd59VegMziwC9CGbcPRf4XVJSikgoTs85nLyCIh579RMeeG4+WZlRju3WPOxYUkWUdYvvd2b2PMGoEW+aWSHwKbAeiACHAZ1jn58ABmgadpHU88MTjiCvoIhnpi7mvj+/T1ZGlD6dm4QdS6qAMt+DcvdV7v5ToClwOTAVWAWsIOh6fiHQ0N1/ouIkkrouOqUD5514BEXFJYx7ci4LlqwNO5JUAYlO+b4deD32IyJVTCQSYeAZncgrKOL16csY+/gcRg3uT+e2DcKOJilMo0KKSEIikQhDzunC93q3Ir+giFGPzGLpyo1hx5IUpgIlIglLS4vw4wu7cWy35uzML2TkxJks/3pL2LEkRalAich+iaZFuOGSHvTp1IStO3Zxx/gZfLVuW9ixJAWpQInIfkuPpnHzFb3o1uEwNm3LZ/ifpvPNhh1hx5IUk3CBMrO+Zva8mS2IDSR7k5ldkMxwIlJxZWZEuf3KPnQ6vD7rN+dx+5+m8+3mnWHHkhSS6HxQPwD+CWwADMgg6AH4tJldlbx4IlKRZWelM/KafhzRsi7fbNjB8PEz2LQ1P+xYkiISbUGNBn7u7kOBQgB3vxv4McGLvCJSRVXPzmD0kP60aVqbVWu3MWLiDLbuKAg7lqSARAtUR+Afe1j+NtD60MURkcqoVvVMRg/tT/PDarBs9RbufHgmO/J2hR1LKrlEC9QqgjH3SjuFYFQJEani6tXKZmzuABrVr86SLzcx+tHZ5BUUhh1LKrFEC9RdwAQzuxWIAmea2W8JBoe9N1nhRKRyaVi3Gnfl5tCgTjaffPEtdz8+h12FRWHHkkoqoQLl7pOAi4HTCCYoHAX0Ay5x90eTF09EKpsmDWowZmgOdWpmMn/JOn41aR6FRcVhx5JKKKGx+ADcfSrBYLEiImVq2bgWY4bmcNtD05n9yRp+9+wH3HBJT6JpkbCjSSWSUIEys+vKWu/uDx2aOCKSKg5vVodRQ/ozfPwM/j3/K7Iyovz4gm6kqUhJghJtQf1yD/s1IuhyPh1IqECZ2SnAPUB7YC1wr7tPMLMmwJ+AE4ASgqk8fuLum2P7XQjcTTDtx7vAle6u8f5FKrgOreox4uq+jHx4Fm/N+ZKszChDzulCJKIiJfuW6DOow0v9tAQaAK+w5+7n/8PMWgIvAGOBugTPtMaZ2WnAI8BGoBnQhqAQjYnt1wl4FLgyds6lwOQEr09EQnZUu4bcPqgP6dE0Xpu2jKfe+DTsSFJJHPBYfO6+DRgJ3JDgLm2AZ9z9JXcvdve5wDvAAILp4nPdfSdQB6hBMHMvBNPJv+ru09w9D7gVGGBm7Q80u4iUrx7WiJuv6EVaWoTn317KX/7hYUeSSuBgB4vtCGQnsqG7v+fuubu/m1l94FhgvrsXuHuBmT1L8F5Vbf5727ATsCjuODuAlUCXg8wuIuWo31FNufGSHkQi8PQbi5ny78/DjiQVXKKdJJ7bw+I6wInAE/t7UjOrQ3B7cDYwJW7VIOBa4DHgReA4oCZQepjkHUD1/T2viITruO4tyC8o4g/PLeCRKR+TnRnltH5two4lFVSinSS2l/peQjBw7HPAU/tzQjPrQFCUFgGXuvt/XpCI3cLLM7ObgKWxVtZ2oFqpw1QHNAGNSCV0St/W5O8qYsJLH/HgXxeSmRHlxJ4tw44lFVBCBcrdBx2Kk5nZcQTFaTxwm7uXxJbPBW5y93/FNs0i6CG4naCQWdwxqgOtiLvtJyKVy5nHtCWvoIgnX1/E/ZPnk5URJadrs7BjSQWT6C2+EYke0N1H7+UY7Qi6j9/u7g+UWv0+MMrMFgIR4DfAJHfPN7NngGlmdgIwExhH8NxqSaKZRKTiOf+k9uTlF/KXfyzh3qfncfugvvTq2DjsWFKBJHqL72jgB8BOYAGQD3QFmsS+775NV0IwNceeDANqEXQtHxe3/EHgRuA+YDFBy+l5gt56uPtHsTmnxgPNCZ5baaJEkRRw6fePJK+giCn//pxxT8zhzsH96XJEw7BjSQURKSkp2edGsYFhWwMDY93LMbNMgp522939Z0lNeRDMrA2w7O2336ZFixZhxxGRUkpKSnjwrwuZOmsF2ZlRxuTmcGTr+mHHknKyatUqTj75ZIDD3X15/LpEu5lfBQzfXZwA3L0A+HVsnYjIAYlEIlx33tGc0LMFeQVF3DlxJp+v2hR2LKkAEi1Q29jze0cDgG8PXRwRqYrS0iL8/Efd6d+lKdvzChkxcSZfrtkSdiwJWaLPoP4APGpmvYAPYsv6A4MJni2JiByUaDSNX17Wi7sen837i9dyx4QZjBt2DM0a1gw7moQk0bH4fg1cT/Bi7qPAAwSdJM5x98eTF09EqpKM9DRuvbIPXY9oyIYt+QwfP4O1G0u/py9VRaLdzG8j6Pb9SJLziEgVl5URZfhVfRkxYQaLV2zkjvEzuGfYMdSrndCoapJCEn0GdTOQmcwgIiK7VctKZ+Tg/rRtXofV67czfMIMNm/LDzuWlLNEC9QrwA1mdlgyw4iI7FazWgajh/SnZeNafLlmKyMfnsn2nbvCjiXlKNEC1Qm4DlhjZlvNbG38TxLziUgVVqdmFmNzc2jaoAafr9rMqEdmsTO/MOxYUk72pxefiEi5q187m7G5Odz84DQ+Xb6BsY/NZuQ1/cjMiIYdTZIs0cFinwQws2oE07WnAZ+7+9YkZhMRAaBR/erclZvDLQ9O48PP1jPuybncdmUfMtIPdko7qcgS+l/XzKJm9iuCadnnE7wLtc7MxptZoq0wEZED1uywmozJzaFW9UzmffoN9/35fYqKive9o1Raif7nx10EU69fAbSM/VxBMIBswiOdi4gcjNZNajN6aH9qZKcz/cPV/OG5BRQX73s8UamcEi1QVwCD3f05d1/t7l+5+3PAUIJZcEVEysURLeoy8pr+ZGdG+ee8lYx/8UMSGfRaKp9EC1RN4LM9LP8C0Nj4IlKuOh5en+FX9SUjPY03Zi7nsVc/UZFKQYkWqLnsecy9HxNMNigiUq6Obn8Yt13Zh/RohJff/ZxnpnrYkeQQS7SDw83AO7FZbWfFlvUD2gDfP/SxRET2rVfHxvzisl78etJcJr/lZGdGOe+k9mHHkkMk0cFi5wHdgbcIOkg0BF4FjnT3mcmLJyJStgFdm/Gzi3oQicATry/i9WlfhB1JDpFEB4sdAfzG3X9RanltM/utu9+QlHQiIgk4qVdL8ncV8dBfFzL+pY/IyozyvT6tw44lB2mvBcrMmgN1Yl9HAv80sw2lNusG5AIqUCISqtP7tyG/oIhHX/mYB55bQFZGOsd2bx52LDkIZbWgegMvAru7xvx7L9s9ekgTiYgcoHOOb0deQSF//vti7nvmfbIyo/Tp3CTsWHKA9voMyt1fJugE0Q6IAH2Aw+N+2gAN3X1w0lOKiCToR9/rwHknHkFRcQnjnpzLfNd41pVVmc+g3P3L2EcNeCUilUIkEmHgGZ3ILyjitenLGPv4HEYP6U/ntg3Cjib7SYVHRFJOJBJh8Dld+F7vVhTsKmLUI7NY8uXGsGPJflKBEpGUlJYW4ccXduO4bs3ZmV/IyIkzWbZ6c9ixZD+oQIlIyoqmRbj+kh707dyEbTt3MWLCTFZ+o1mCKov9LlBmVs/MVNhEpFJIj6Zx0+W96NbhMDZty+eOCTNY8+32sGNJAhIuNGZ2k5mtA9YBbczsSTP7o5llJC+eiMjBy8yIcvugPnRu24BvN+cxfPwM1m/aGXYs2YdEJyz8JXAt8FMgP7b4r8APgXHJiSYicuhkZ6Yz4uq+tG9Zl2827GD4+Bls3JoXdiwpQ6ItqGuAXHd/FigGcPdXgYHAxUnKJiJySFXPzmDUkP60aVqbr9ZtY8SEmWzdURB2LNmLRAtUK2DpHpZ/CdQ7dHFERJKrVvVMxgzNoflhNVn+9RZGTpzJjrxdYceSPUi0QL0PXBT3fffwR8OADw5pIhGRJKtbK4u7rs2hcf3qLF25idGPziYvvzDsWFJKogXqRuCXZjYVyALuMrMPgKuAm5IVTkQkWRrUqcbY3Bwa1Mnmky++5a4n5rCrsCjsWBInoek23H22mRlBi2kLUA34O3CWu3+V6MnM7BTgHqA9sBa4190nmFkj4PfAyQTj/r0B/MzdN8b2mwRcCMT/J05Xd9fELyJywJo0qMHY3BxufXA6C5as41eT5nHLwN6kR/UmTUWQ6Iy6uPtagmk3DoiZtQReIOhYMQXoCUw1s+UEhW8zwSC0GcBTwIPAJbHdewDnuPvfD/T8IiJ70qJRLUYP7c9tD01n9idr+O0zH3DjpT2JpkXCjlbllTUf1HOJHsTdL0xgszbAM+7+Uuz7XDN7BxhA0DNwlLtvj537YeCPsc/VgCOBBYnmERHZH4c3q8OoIf0ZPn4G7y34iqyMKD+5sBtpKlKhKqsFdUhftXb394D3dn83s/rAscBT7j6i1ObnAPNjn7sR3Np72Mz6ASuBEe7+2qHMJyJVW4dW9Rh5TT9GTJzJP+Z+SXZmlCE/7EIkoiIVlr0WKHcflKyTmlkd4BVgNsHtvvh1vyAoUDmxRbUICtsoYCFwNvCcmfV394XJyigiVU/ntg0YPqgPox+dzWvTl5GVGWXgGZ1UpEKS8DMoM8shGE2iM1AALCLo5PDp/pzQzDoQFKVFwKXuXhxbngE8AJwFnOTuiwHc/U3gzbhDvGBmgwgKlQqUiBxS3a0Rt1zRi3FPzuWFf31Gtax0fnSKhR2rSkp0qKMLCKZ8rwE8D7wKHAYsMLOzEj2ZmR1H0Gp6GTjf3fNiy2sBbxFMM9/H3RfE7XOWmQ0sdahMQGOUiEhS9D2qKTde0pO0CDz998W8/O7nYUeqkhJtQd0DXO/uD8QvjN2Ou5egYJXJzNoBrwG3lz4OMJmgWB7r7jtKrYsCvzezTwleGP4Rwe2/axLMLiKy347t3pz8XYX8/i8LePSVj8nOjPL9/m3CjlWlJFqgmgJT97D8VWBMgscYRvA8aZyZxQ8w+xbwA4JBaNcGr1sBsMndW7j7y2Z2O/As0ARYDJwZNx29iEhSfK9Pa/ILihj/0kc89MJCsjKjnNizZdixqoxEC9SLBM+fri+1fCDwt0QO4O43ADckHu07+z5I8F6UiEi5OuOYtuQVFPHE64u4/9kPyMyIMqBrs7BjVQmJvgdVA7jEzE4jeIZUBHQFugNPJzWhiEjIzjupPTsLCvnLW0v4zdPzyBrUl14dG4cdK+WV1Ulie9zPWuBJguIEwXOhT1BxEpEq4tLTjuSc49tRWFTCuCfm8OFn68KOlPJCeQ9KRKSyiUQiXHVWZ/ILinhj5nLGPDqb0UNy6Hh4/bCjpaz9eQ+qO9CJoPUEwaCuWUBPdx+ahGwiIhVKJBIh99yu5BUU8q/3VzHqkZmMvXYAR7SoG3a0lJRQgYr1ohsDbCN4HrUZqBNbnVAnCRGRVJCWFuFnP+pOwa5ipn+4mhETZjJu2ABaN6kddrSUk+iY8kOBX7p7beBrgg4SzYFZwNwkZRMRqZCi0TRuvLQnvTo2ZuuOAu4YP4PV67aFHSvlJFqgmhBMlQHBqOL93X0NwWSFlycjmIhIRZaRnsYtA3vT9YiGbNyaz+3jZ7B2Q+lxBuRgJFqg1gENYp+XAEfHPn8F6IUAEamSsjKiDL+qLx3b1Gf9pp0MnzCDDVs0CtuhkmiBmgJMNLNuwL+AK8zseIIXb1ckK5yISEVXLSudEdf0o12LOny9fjvDx89g87b8sGOlhEQL1C8InjUdRTCe3j+BtwlGkvhFcqKJiFQONatlMGpwf1o1qcXKb7YyYuJMtu3cFXasSi+hAuXuO9w9192fdvcSd78SqAs0cHf14hORKq9OzSzGDM2hacMafPHVZkY9PJOd+YVhx6rUyhrq6AeJHMDMUJESEYH6tbMZm5vDLQ9OY/GKjYx9bDYjrulHVkZ03zvL/yjrPahEp1Qv4b8v74qIVGmN6lVnbG4Otz44jQ8/W889T87ltiv7kJGe6BMV2a2soY70tykicgCaNazJmKE53PrQdOZ9+g2/+fM8brqsF9Gofq3uj7Ju8XUCFrt7cezz3pTs77TvIiKprlWT2owe0p/bx89gxodfc/9f5nP9RT1IS4uEHa3SKKucfww0jPv8UfXgM+wAABWMSURBVOzPPf2IiEgp7VrU5c7B/cjOjPLO+6v404sfUlJSEnasSqOsAnU4wQu6uz+3jf1Z+qdtMgOKiFRmR7auz4ir+5GZnsbfZy7n0Vc+UZFKUFnPoFbs6XNpZqaRJEREytDliIbcNqgPYx+bzZR/f052VpTLvt8x7FgVXqKjmXcA7mXP0200SvQ4IiJVVc8jG3PT5b24Z9I8/vLWErIz0zn/pPZhx6rQEu1S8iegVezPZsBDwBsEz6gGJyeaiEhq6d+lGddf1J1IBJ58fRGvvvdF2JEqtEQLVD9gqLv/FlgIzHL364CbgUuSFU5EJNWc0LMlw84Pxtue+PJHvDVbw5nuTaIFKgKsiX1eDHSPfX4Z6HGoQ4mIpLLT+rVh8P8dBcADzy/g3/NXhZyoYkq0QC0Ezo19/gQ4Pva5OUHxEhGR/XD2ce247PQjKSmB+575gFkffx12pAon0QJ1J3CPmQ0DngZONbN3gL+iKd9FRA7Ij75nXHBye4qLS/jVpHl84GvDjlShJDqa+VSgA/CGu68GcoA5wO9QJwkRkQN2+ekdOevYthQWFXPX43P4+PP1YUeqMBLtZn4L8Iy7fwHg7h8TTPcuIiIHIRKJcM3ZR5GXX8hbc75k9KOzGJs7gA6t6oUdLXSJ3uK7CPjCzN4zs6FmVj+ZoUREqpK0tAjDLujG8d1bsDO/iBETZ7Js9eawY4Uu0Vt83YDOwFvAT4GvzexVM7vYzKolM6CISFUQTYvw84u70++oJmzfuYs7Jsxg5Tdbw44VqoTHfvfAaHfvDPQCFgATgG+SFU5EpCpJj6Zx0+W96GGN2LytgOHjZ7Dm2+1hxwrNfk1OYmaZZnY2wQu6w4CNBKNKiIjIIZCRHuXWK3vTuW0DNmzJ4/bxM1i/aWfYsUKRUIEyszPNbBKwFngM2Aac4+6t3f2WZAYUEalqsjPTGXF1Xzq0qsvaDTsYPn46G7fmhR2r3CXagvoLkAlcDjRx91x3/7eZ1TOzHycvnohI1VQ9O4NRg/tzeLPafLVuOyMmzGTL9oKwY5WrRAtUY3e/yN1fdfdCMzvVzCYDq4H7k5hPRKTKqlk9k9FDcmjRqCbLv97CyIdnsiNvV9ixyk1C70G5+zYzawMMAq4EWhDc5nsY+GOiJzOzU4B7gPYEtwvvdfcJZtYI+D1wMsHQSW8AP3P3jbH9LgTuBpoC7wJXuuuVaxFJfXVrZTE2N4dbHpzGZys3MeqRWYwa3J/srNSf5ajMFpSZZZnZpWb2NvAZcCvgQAlwnLv/1N2XJHIiM2sJvACMBeoCFwPjzOw04BGgkGCG3vZAPeDB2H6dgEcJCmMDYCkwef8uU0Sk8mpQpxpjcwfQsE42i5Zt4K7H51CwqyjsWEm31wJlZg8RjGA+EdgKXE1wq+9UggK1v+3MNgSjUbzk7sXuPhd4BxgAFAOj3H27u28iaJkdE9vvMuBVd5/m7nkERXKAmWmmLxGpMhrXr87YawdQt1YWC5au41eT5lFYVBx2rKQqqwWVS/COUy5wtbs/ufuW24Fw9/fcPXf399hoFMcC8939HHf/LG7zc4D5sc+dgEVxx9kBrAS6HGgWEZHKqPlhNRkzNIda1TOYs2gN9/35fYqKS8KOlTRlFagTgX8TPBtaY2bvmNlPzKzFwZ7UzOoArwCzgSml1v2CoEDdHFtUE9hR6hA7gOoHm0NEpLJp07Q2o4b0p3p2OtMWruaB5+ZTnKJFaq8Fyt3fdfchQBOCsfg2APcCK2L7nWtmtff3hGbWAZhF0Do7392LY8szzGw8cD1wkrsvju2yHSg9nFJ1gk4aIiJVTvuW9RhxdT+yMqO8PXclE176kJKS1CtS++xm7u4F7v6Cu59LUKyuA6YDo4DVZvZIoiczs+MIWk0vExSnvNjyWgTj/PUG+rj7grjdFgEWd4zqQCvibvuJiFQ1nds2YPigPmSkp/G3Gct54rVFKVek9qufYqwDwwRggpm1Inhx95JE9jWzdsBrwO3u/kCp1ZMJiuWxsWdM8Z4BppnZCcBMYBzBc6uEeg+KiKSqbh0accsVvbn7iTm8+M5nZGelc/Gptu8dK4kD7kjv7l8Cd8V+EjEMqEXQtXxc3PK3gB8A+cBas//85W5y9xbu/pGZXQWMJ5hifjZwwYHmFhFJJX06N+HGS3vym6fn8czUxWRnRvnhCUeEHeuQKLc3vdz9BuCGA9z3BYJ3qEREpJRjuzWnYFcR90+ez2OvfkJ2ZpTTcw4PO9ZB26/RzEVEpGI6uXcrcs/tCsBDL3zIP+d9GXKig6cCJSKSIs4YcDiDzuwEwO8nz2f6wtUhJzo4KlAiIink3BPbc/GpRnEJ3Pv0POYuWhN2pAOmAiUikmIuPtX44QlHUFRcwrgn57JwybqwIx0QFSgRkRQTiUQYdGYnTs9pw67CYsY8PptFy74NO9Z+U4ESEUlBkUiE3B925aReLckvKGLUI7P4bOWmsGPtFxUoEZEUlZYW4acXdmPA0c3YkVfIiIkzWPH1lrBjJUwFSkQkhUWjadx4SU96dWzM1h27GD5hBqvXVY6hTFWgRERSXEZ6GrcO7E3XIxqyaWs+t4+fwdoNpUeVq3hUoEREqoDMjCjDr+pLxzb1Wb9pJ8PHz+DbzTvDjlUmFSgRkSqiWlY6I6/pxxEt6vD1t9u5Y8IMNm/LDzvWXqlAiYhUITWqZTBqSA6tm9Ri5TfbGDFhJtt2FIQda49UoEREqpjaNTIZMzSHZg1r8MXqzdz5yCx25O0KO9b/UIESEamC6tXOZmzuABrVq4av2MjYx+aQv6so7FjfoQIlIlJFHVavGmNzB1C/dhYffb6eu5+Yw67CilOkVKBERKqwpg1rMGZoDrVrZPLB4rXc+/T7FBUVhx0LUIESEanyWjWpzZihOdSolsHMj77md8/Op6i4JOxYKlAiIgJtm9dh1OB+VMuK8u78VfzphYWUlIRbpFSgREQEAGtdnzuu7kdmehpTZ63gkSkfh1qkVKBEROQ/urRryG2D+pAejfDKe1/w9N8Xh5ZFBUpERL6j55GNueny3qSlRXjuH0t4/u0loeRQgRIRkf/Rv0tTrr+4B5EITPrbp7zy78/LPYMKlIiI7NEJPVow7PxuADw85WOmzlpRrudXgRIRkb06rV9rBp9zFAAP/nUB73ywqtzOrQIlIiJlOvvYdlzxg46UlMDvnv2AmR+tLpfzqkCJiMg+XXByBy44uT3FxSX8+ql5vL/4m6SfUwVKREQScvnpHTn72LYUFpVw9+Nz+Ojz9Uk9nwqUiIgkJBKJcM3/HcWpfVtTUFjMmEdnsXjFhqSdTwVKREQSFolEuO78ozm+ewt25hcx6uHkzSWVnpSjiohIyoqmRbj+4u5kpKexeMUGIpFIUs6jAiUiIvstGk3jZxd1T+o5dItPREQqJBUoERGpkMr1Fp+ZnQLcA7QH1gL3uvuEuPXZwDvAPe7+ctzyScCFQGHc4bq6+xflkVtERMpfuRUoM2sJvAAMBKYAPYGpZrbc3aeaWVdgItB3D7v3AM5x97+XV14REQlXebag2gDPuPtLse9zzewdYICZLQPeBu4GmsbvZGbVgCOBBeUXVUREwlZuz6Dc/T13z9393czqA8cC84HVQFt3/x1QevrGbgS39h42s3Vm9oGZnVleuUVEJByhdJIwszrAK8BsYIq7b3P3rXvZvBbwHjAKaAbcBTxnZkeXS1gREQlFub8HZWYdCJ5BLQIudffisrZ39zeBN+MWvWBmg4CzgYUJnDIKsGbNmgMLLCIiSRP3uzlael159+I7jqA4jQduc/fSt/P2tM9ZQH13fzJucSaQl+BpmwJceuml+5lWRETKUVPgO9P2lmcvvnbAa8Dt7v7AfuwaBX5vZp8C7wM/AnKAaxLcfy7Bs66vgaL9OK+IiCRflKA4zS29ojxbUMMInieNM7NxccsfdPeb97aTu79sZrcDzwJNgMXAme7+ZSIndfd8YNqBxxYRkST7fE8LIyUl+7zLJiIiUu401JGIiFRIKlAiIlIhqUCJiEiFpAIlIiIVkgqUiIhUSCpQIiJSIalAiYhIhVTuY/FVJrEBaccDXYEvgKvc/X/edq7MzKwP8Jq7N4p9zwT+CJxPMPLGb919XNz2F/LfaVHeBa5097XlHvwA7W3SzFS/boDYLAB3A4cTXPuvq8q1A5hZXeBDYIS7P1EVrtvMrgImAPlxi4cRDHxQ4a9dLai9iP3jnQL8BahLMIr6m2ZWO9Rgh4iZRczsGoKBeDPjVo0CDGgH9AYGmtkVsX06AY8CVwINgKXA5HKMfVDiJs0cS/C/6cUEI5ucRgpfN4CZNQX+Ctzs7rWAC4D7zawHKX7tccYDzeO+V4Xr7gHc5+41436epJJcuwrU3p0AZLj7/e6+y90nA58QjAWYCkYB1xL8so43ELjL3Te6+3LgN8DQ2LrLgFfdfZq75wG3Ekw42b6cMh+sNsQmzXT34lhr+B1gAKl93bj718Bh7v6GmaUR/OIpBLaS4tcOYGYDgdrAR3GLU/66CWYu39Nkr5Xi2lWg9q4T8GmpZYuBLiFkSYbx7t4TmLd7QewWSFOCqVB2i7/mTvHr3H0HsJJK8neyj0kzU/a6d3P3rWZWneB2z5vAg8A6UvzazexwYCRwVdyylP63DmBmUYLHE5eb2Woz+8zMbjGzelSSa9czqL2rCewotWwHUD2ELIecu6/ew+KasT/jrzv+mlPm76TUpJnvxxan/HUTTFNTg+AX19+AnbHlKXntsV/STwO/cPc1ZrZ7VVX4t34YwX+APgmcC3QkeGyx+5Z+hb92Fai92w5UK7WsOrAthCzlZXvsz/jrjr/mlPg7KT1pJv+9ppS+boDYBKEFwDwzmwj0iq1K1Wu/A3B3f7HU8pT/t+7ua4Dj4xYtMLMHgNNj3yv8tesW394tIniIGO9IvtssTinuvhFYw3evO/6av/N3Ertd1IpK9HcSmzRzNvAycL6751WR6z7ezN4vtTgLSPVrvwg438w2mdkmgttUDxF0ekrl68bMOpvZqFKLd0/2WimuXS2ovfsXEDGz6wm6Y55HcFvkpVBTJd9TwEgz+5Cgqf8L4Pexdc8A08zsBGAmMA6Y7+5Lwgi6v/YxaWbKXnfMAqC5md1AcF19gauBHxL8skrJa3f3I+O/m9kC4P5YN/NtpOh1x2wCbjSzVQS98roDPwV+TNDhq8Jfu1pQe+HuBQRN4fOADcDtwDnuvi7UYMk3AviY4B/wXIJu2eMB3P0jggfN44H1QGeC7sqVRfykmdvifn5Fal837r4Z+AHBs4gNwETgGnd/lxS/9jKk9HW7+1fA2QS987YQXN8Yd/8rleTaNWGhiIhUSGpBiYhIhaQCJSIiFZIKlIiIVEgqUCIiUiGpQImISIWkAiUiIhWSXtSVlGFmy4HWBCNEvFBqXX3gG2CzuzdM8HhdgQbu/q9DHLVCiM0P9aq7R8rYphYwHTjB3TeUQ6YWwKtAP3fP39f2ktrUgpJUsws4Zw/Lzwai+3msKQQvKVZlY4Fny6M4Abj7KuBt4JbyOJ9UbCpQkmreAc4ws9J3B84lGLZlf+y1ZVEVmFkTguGQxpfzqf8A3JAqk4PKgdMtPkk1bwI5wHHAPwHMrCZwEsEkjTfv3jA2y+yvgdOAbODvwM/c/Wsze4fgduEDZna+u5+wj+3bAMsIRs++Hpjt7j+IDxYrmr8jGDamDsE0Hze4+5y4/S8hGMj0MGAqMNTdv43t34Tgl/fpBCNLvw7cGBvGKJH1RxAUmxxgCcGYa2XJBWbEBtPFzK4kGMttEUGL9NcEv0N6xLJfBWwGfuPu98ft82OCKS9uIxig9k/AJILx4Y4mGCfwUndfBuDuX5rZUmAQ/x0fTqogtaAk1eQBbwD/F7fsDILxxv4zjqKZZRDcSmpFMEbdSQTTgb9sZhGCFtcqgl+q5yaw/W5nAv2BX+4h20+A7wNnEdw6XAL8tdT+4whmOj4BOBx4Pm7diwT/n+0fO0Y7vjsV917Xx/L/jWAqhV7AaPZ9G+0MgiIcrzvBLLw9CIoMBAW7FdCPYPbVcbEZbHfrAhxDMDnkzbGfVwhuHw4A6hP8x0O8v/PfaSGkilILSlLRiwS/6H8W+34uwWCY8U4DjgC+t3vyRjP7EUFL4Hvu/paZFQFb3X1DrEPBXrcHlsaOe38Zoz63IZj4bbm7r4uNLN6d7/6H4h3uPjV2/KuB+WbWEWhCMJr+ibs7D5jZpcBXZnYUQYurrPUtiRWR2POkRbF5scbtKWhsor/uwPA9rB4Vm2uI2ASAO4Er3H0r8ImZ9QSuI5goD4IpHq5z97WAm9l9wGR3fz12jMkEBTXeIoIWnFRhakFJKnodaGpm3cwsm6DVUnqalM7AiviZhWMP6Jez544RiW7/eRm5/khwa2+1mb1L8HznQ3cvitvmvbjPCwmmZz8qdo7qwLe7R2EnaIFBMHfPvtYfFcsf39lhThlZGxB0KllfavmO3cUpzoJYcYo/7lFx37fEitNuO4Ev4r7nEdz6i/ctUD9WKKWKUoGSlOPuWwhux50DnAJ8Ept6IN7O/9kxEGHP/79IdPu9bYe7LyUoFhcADtxI0EJqErdZ4R6OX0Rwt2MF0K3UT3uC5277Wr/7WPEK9pYVKN7LPnu6vtKZ02KZd9tVxvH3JhrbRtMtVGEqUJKqXiQoUHu6vQfwKdA61vEBADNrRtAxYnFsUcl+bl8mMxtMMKfYy+4+BOgANCbo0LFbj7jP3Qlujy2Mnb8ZwS3Hz9z9M4ICcx/QKIH1HwJtzKzxXs5V2rcEheWwBC7tKDPLjPveJ5b5YDQE1sWmqJcqSs+gJFVNIeix1pagQ0Bp/yD4JTo59iwI4LcEt8X+Efu+DehoZo0S2L5ZApnqAmPMbAPBM6szCFoKC+K2uS+2Ph+YALzm7p+b2TKCyeUmm9kvCFotDxB0MFhO8CysrPUrCIrYpNj6lgQdGvbI3UvMbD5BL7s39nFdjYGHzOxego4Sg4HLEvj7KMvRBB1bpApTC0pSUmzm4+nA0t3dl0utLyFoYa0jeHfqbWA1cHJsNmUIujhfDkxNcPt9+S1Bz7fHCW7xDQXOK9Wp4nHgWf5bEC+J5S0m6Nq9AfgX8C5BK+cH7l6UwPpCgt6Hu4BZBN3d79tH3teB4xO4rg8JCuoHBDO1Xld6JI8DcBzw2kEeQyo5zagrUgHEvQfVxd0/DjkOAGbWnKCQti3VySF+mzuBM9291yE8bwdgNtCqVOcLqWLUghKRPYp1LHmSoKVXnoYBf1BxEhUoESnLbcBFZtagPE4WGyz2RPbyfpZULbrFJyIiFZJaUCIiUiGpQImISIWkAiUiIhWSCpSIiFRIKlAiIlIh/T8NCUZ4iERw1gAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "xs = linspace(0, 525, 21) * radian / s\n",
    "taus = [compute_torque(x, system) for x in xs]\n",
    "plot(xs, taus)\n",
    "decorate(xlabel='Motor speed (rpm)',\n",
    "         ylabel='Available torque (N m)')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Simulation\n",
    "\n",
    "Here's the slope function that computes the maximum possible acceleration of the car as a function of it current speed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "def slope_func(state, t, system):\n",
    "    \"\"\"Computes the derivatives of the state variables.\n",
    "    \n",
    "    state: State object\n",
    "    t: time\n",
    "    system: System object \n",
    "    \n",
    "    returns: sequence of derivatives\n",
    "    \"\"\"\n",
    "    x, v = state\n",
    "    r_wheel, speed_ratio = system.r_wheel, system.speed_ratio\n",
    "    mass = system.mass\n",
    "    \n",
    "    # use velocity, v, to compute angular velocity of the wheel\n",
    "    omega2 = v / r_wheel\n",
    "    \n",
    "    # use the speed ratio to compute motor speed\n",
    "    omega1 = omega2 / speed_ratio\n",
    "    \n",
    "    # look up motor speed to get maximum torque at the motor\n",
    "    tau1 = compute_torque(omega1, system)\n",
    "    \n",
    "    # compute the corresponding torque at the axle\n",
    "    tau2 = tau1 / speed_ratio\n",
    "    \n",
    "    # compute the force of the wheel on the ground\n",
    "    F = tau2 / r_wheel\n",
    "    \n",
    "    # compute acceleration\n",
    "    a = F/mass\n",
    "\n",
    "    return v, a      "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Testing `slope_func` at linear velocity 10 m/s."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>values</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>x</th>\n",
       "      <td>0 meter</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>v</th>\n",
       "      <td>10.0 meter / second</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "x                0 meter\n",
       "v    10.0 meter / second\n",
       "dtype: object"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "test_state = State(x=0*m, v=10*m/s)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(10.0 <Unit('meter / second')>, 14.201183431952662 <Unit('newton / kilogram')>)"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "slope_func(test_state, 0*s, system)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can run the simulation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>values</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>success</th>\n",
       "      <td>True</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>message</th>\n",
       "      <td>The solver successfully reached the end of the...</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "success                                                 True\n",
       "message    The solver successfully reached the end of the...\n",
       "dtype: object"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "results, details = run_ode_solver(system, slope_func)\n",
    "details"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And look at the results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>x</th>\n",
       "      <th>v</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>2.88</th>\n",
       "      <td>57.32626824256294 meter</td>\n",
       "      <td>38.632749274430566 meter / second</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2.91</th>\n",
       "      <td>58.49067227876934 meter</td>\n",
       "      <td>38.99375142852514 meter / second</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2.94</th>\n",
       "      <td>59.66589334398368 meter</td>\n",
       "      <td>39.353885587648904 meter / second</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2.97</th>\n",
       "      <td>60.851905429699734 meter</td>\n",
       "      <td>39.713153838812325 meter / second</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3.00</th>\n",
       "      <td>62.0486825899462 meter</td>\n",
       "      <td>40.071558264007834 meter / second</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                             x                                  v\n",
       "2.88   57.32626824256294 meter  38.632749274430566 meter / second\n",
       "2.91   58.49067227876934 meter   38.99375142852514 meter / second\n",
       "2.94   59.66589334398368 meter  39.353885587648904 meter / second\n",
       "2.97  60.851905429699734 meter  39.713153838812325 meter / second\n",
       "3.00    62.0486825899462 meter  40.071558264007834 meter / second"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "results.tail()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "After 3 seconds, the vehicle could be at 40 meters per second, in theory, which is 144 kph."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "40.071558264007834 meter/second"
      ],
      "text/latex": [
       "$40.071558264007834\\ \\frac{\\mathrm{meter}}{\\mathrm{second}}$"
      ],
      "text/plain": [
       "40.071558264007834 <Unit('meter / second')>"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "v_final = get_last_value(results.v)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "144.2576097504282 kilometer/hour"
      ],
      "text/latex": [
       "$144.2576097504282\\ \\frac{\\mathrm{kilometer}}{\\mathrm{hour}}$"
      ],
      "text/plain": [
       "144.2576097504282 <Unit('kilometer / hour')>"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "v_final.to(km/hour)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plotting `x`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO3deXxU1d3H8U9WkhAS1kACyCYcRBbZQSogVluXalXqrqhQ0ccu2lq7aEHcW5dqtX0Qd23VukvVuoNSWUT29bDvJBCW7Mtk5j5/3IEnRAgTSObOTL7v1yuvZO6cmfu7Xsk3586558Q5joOIiEikife6ABERkcNRQImISERSQImISERSQImISERK9LqA42WMaQIMBnYCfo/LERGRukkAsoH51tqK6k9EfUDhhtMsr4sQEZHjchrw3+obYiGgdgL885//pF27dl7XIiIidZCbm8uVV14Jwd/l1cVCQPkB2rVrR4cOHbyuRUREjs13PqLRIAkREYlICigREYlICigREYlICigREYlIsTBI4qgKCwvZtWsXPp/P61LqRdOmTenQoQPx8fr7QkRiV8wHVGFhIXl5ebRv357U1FTi4uK8Lum4BAIBtm/fTn5+PllZWV6XIyLSYGL+T/Bdu3bRvn170tLSoj6cAOLj42nbti0FBQVelyIijVzunhKKSisb7P1jPqB8Ph+pqalel1GvkpKSqKqq8roMEWnEFqzOY+KDn/Pgi/MbbB8xH1BATPScqou14xGR6LIlt5A/v/wtgYBD3+6tG2w/jSKgRESkfhQUV3D3s/MoLa9iRL8cfjKmR4PtSwElIiIh8VX5ue/5b8jbW0r3js255bL+xMc33BWdsI7iM8ZkA/8LnA6UA9OstX80xiQDTwJjcedjetRa+0A4axMRkSNzHIcnXl/Mqk17aZ2Zwp3XDyUluWEjJNw9qPdwZ6xtCwwDxhljrgCmAAbohrt8xjhjzDVhrs1Tb7/9NgMHDiQvLw9wZ2cfPXq0RuuJSER4/bM1zFiwjZTkBO68figtM1IafJ9h60EZY4YCXYER1lofsNEYMxooAx4GrrXW7gP2GWMeBiYCL9V3HVOemcu3q/Lq+22PaNBJbZk8YdhR21100UV8/vnnTJkyhd/85jc88sgjTJ06lczMzDBUKSJyZLMWb+cfH60mLg5+c9UgunVoHpb9hrMHNRBYBtxljNlujFkPXIgbUNnAymptVwN9wlhbRLj77rtZtGgR1113HVdffTVDhgzxuiQRaeTs5r089upCAK7/0ckMOTl86+6F8zOolrgrJn6J25PqCXwE7A4+X1qtbSmQ1hBFhNKb8UqrVq0466yzeO211zj//PO9LkdEGrm8vaXc+9w3VFYF+MGwTlwwsltY9x/OHlQFUGitvctaW2GtXQI8A4wLPl/9bto0oDiMtUWEZcuWMX36dM4991zuvPNOAoGA1yWJSCNVXOZjyjNz2V9cQb/urbnxor5hvwcznAG1GkgLjtg7IBHYB+TiDpI4oCeHXvKLeeXl5dx+++3ceOON3H///ezZs4dnnnnG67JEpBGq8gf404vz2ZpXRMe2zfjduCEkJoT/rqRwXuL7FPdy3iPGmF/jBtJ44CZgAzDZGLMUSAduAx4PY22ee+SRR0hJSWH8+PEkJiZyzz33MGHCBEaOHEnPnj29Lk9EGgnHcZj69lIWr91N8/QmTJ4wjPTUJE9qCVtAWWvLjTGjgCdwh5qXA3+21r5ljPkAeARYgdurmwZMDVdtkeCOO+445PHQoUNZtmyZR9WISGP11ox1fDx3M8mJ8dx5/RDatmyQ4QAhCeuNutbaDcC5h9leDtwc/BIREQ/MWrydFz9wP1351RUDMZ1aelqPpjoSERFWbdzLX4LDya8772RG9MvxuCIFlIhIo7cjv5h7npuHryrA2cM7c+Ho8A4nPxIFlIhII1ZQXMFdT8+lqLSSgT2zmHhhn4hZ0qdRBJTjOF6XUK9i7XhExBsVPnd28p35JXTNyeT2qweR4MFw8iOJnEoaSFJSEmVlZV6XUa98Ph+JiWEd3yIiMSYQcPjLKwvd2cmbpzJpwlDSUrwZTn4kMR9QWVlZbN++ndLS0pjoeQQCAfLy8jSJrIgcl+ffX8HXS3eQlpLIXROG0Soz9egvCrOY/zM8IyMDgB07duDz+Tyupn40bdqU1q0bbpllEYlt02et590v15MQH8cfxg2hU3aG1yUdVswHFLghdSCoREQas9lLd/DMe8sB+MWl/enXo43HFR1ZzF/iExER16qNe3nknwtwHLjq7J6MGdTR65JqpYASEWkEtu0q4p7n5h1cOuOSM3p4XdJRKaBERGLcvsJyJgfvdRp0Ultu8mDpjGOhgBIRiWGl5T7uemYuu/aW0r1jc34bYfc61SY6qhQRkTqr8gd48MX5bNheQHarpkwaP4yUJtEzNk4BJSISgxzH4a//WsSiNbvJTE/mrhuG0bxZE6/LqhMFlIhIDHrxg5XMWLCNJskJTBo/jJzW6V6XVGcKKBGRGDP9q/W8NWMdCfFx/H7cYHqc0MLrko6JAkpEJIbMWrSdZ6YfuBH3FAb2bOtxRcdOASUiEiOWrNnNo6+6N+KOO7cXYwad4HVJx0UBJSISA9Zt2899L8yjyu9w/mldufj0E70u6bgpoEREotyO/GKmPD2Xsgo/I/u3Z/z5vaPiRtyjUUCJiESxfYXlTJ42h/3FFZzSvQ23XDaA+PjoDydQQImIRK3iMh+Tn55D7p5STuyQye+vHUxSYuz8Wo+dIxERaUQqfH7ufW4eG3cU0r5NUyZPGB5xK+Ier7DOeWGMuR54Cqiotvlm4FXgSWAs4AcetdY+EM7aRESihd8f4KGXv2XFhj20zEjh7htOjbpZIkIR7kmZBgCPWGt/V32jMeYBwADdgEzgI2PMdmvtS2GuT0QkojmOw5NvLGHeilzSU5O4e+JwslqmeV1Wgwj3Jb6BwOLDbB8H3Get3Wet3QQ8DEwMZ2EiIpHOcRye+/cKPpu/heQkdwqjTu1id7XwsPWgjDEJQF/gamPMo0Ap8AzuJb9sYGW15quBPuGqTUQkGrz5xVre/XI9iQlx/OHawZzUpaXXJTWocF7iawN8C7wIXAScBLwHJAefL63WthSIzT6riMgx+GjOJl76cBVxcfCrywdG9RRGoQpbQFlrc4FR1TYtNsY8AZwdfJxa7bk0oDhctYmIRLJZi7bz97eWAHDTRX05rX97jysKj7B9BmWMOdkYM6XG5mSgHMjFHSRxQE8OveQnItIofbsqj0decefXu+rsnpx9ahevSwqbcF7i2w/82hizDXgW6A/8AvgZsAKYbIxZCqQDtwGPh7E2EZGIs2LDHh54cT7+gMOFo0/kkjN6eF1SWIWtB2Wt3Q6cjzs6rxB4C7jHWvsmMAlYjhtU84PPTQ1XbSIikWb9tv3c/excKn1+zhraievO6xUT8+vVRVjvg7LWfgEMOsz2ctwbdm8OZz0iIpFoa14Rk6bNobS8ihH9cvifsf0aXTiBpjoSEYkouXtKuHPqbApLKhnYM4tfXzGQhBiZ/LWuFFAiIhFiT0EZf3xqNnsLyzm5ayt+Ny62Jn+tq8Z75CIiEaSguII/PjXbnZm8Y3MmjR9KSnK4Z6OLLAooERGPFZf5mPTUHLbmFdOpXTOm/DT2ZiY/FgooEREPlZb7uOvpOWzYUUBO66bcM/FUMpomH/2FjYACSkTEI+6aTt9gN+8jq0Uq9944ghYZKV6XFTEUUCIiHvBV+bn/+W9Ytj6flhlNuPfGEbRpkXr0FzYiCigRkTCr8gf400vfstDuIjM9mXtvHEF266ZelxVxFFAiImHk9wd4+J8LDi44eM/EU+nYtpnXZUUkBZSISJj4Aw6P/WsRXy/ZQVpKIndPHE6XnEyvy4pYCigRkTAIBBz+9sZiZi7YRkpyAndNGE73ji28LiuiKaBERBqY4zj879tL+fQbd6n2yROGxfxquPVBASUi0oAcx2Hau8v4aM4mkhPjmXT9UHp3a+11WVFBASUi0kAcx+GZ6ct5/78bSUyI5w/XDaFfjzZelxU1FFAiIg3AcRye+/cKpn+1gcSEOP5w7WAG9mzrdVlRRQElIlLPHMfhxQ9W8u6X60mIj+N31wxmcK92XpcVdRRQIiL1yHEcXv7PKt6asY6E+Dh+e80ghvbO9rqsqKSAEhGpJwfC6Y3P1xIfH8dvrh7E8D45XpcVtRRQIiL1wHEc/vHR6oPhdPvVgxjRV+F0PBRQIiLHyXEc/vnRal7/bI0bTlcpnOpDyMs1GmPaAgOBLMAP5AILrbV7Gqg2EZGI953LelcNZEQ/hVN9qDWgjDGJwBXALUA/oBLYByQALYNt5gF/B16z1gYatFoRkQhyYLTeWzPW/X/PSeFUb454ic8YMwpYClwDPAv0ANKstTnW2rZAMtAfeAX4GbDaGDO6wSsWEYkAjuPwwvsrD47Wu/1qhVN9q60H9WvgUmvtssM9aa11gOXBr78bY/oDdwMza9uhMaY5bvBNsta+YIxJBp4ExuJeOnzUWvtAXQ9ERCRcHMfhmfeWM33WhoNDyTVar/4dMaCstefX5Y2stYuAH4XQdCrQvtrjKYABugGZwEfGmO3W2pfqsn8RkXAIBNy59T74eiOJCXH89prBDNN9Tg2iLoMk0oAuQJOaz1lrF4b4HuOADKB6r2wccK21dh+wzxjzMDARUECJSEQJBBz+/tYSPp672Z1b71rNENGQQgooY8xVuD2fVCCuxtMO7qCJo71HF2AycCrwUXBbcyAbWFmt6WqgTyh1iYiEiz/g8MTri/h8/laSE+O547qhDOiZ5XVZMS3U+6AewB0o0RU3UKp/HfXCqzEmAfgHcJu1NrfaU+nB76XVtpUCaSHWJSLS4Pz+AI++soDP52+lSXICk8YPUziFQaiX+DKAJ621m49xP38ErLX27RrbS4LfU6ttSwOKj3E/IiL1ylcV4OF/fsvspTtJbZLA5AnDOblrK6/LahRCDaiXgWuBO45xP5cBOcaYi4KPm+HeOzUE94ZfA2wPPteTQy/5iYh4otLn54EX5/PtqjyapiRy1w3D6dlJK+GGS6gB9RCw0BhzJbAJOOSGXGvtmNpebK3tWf2xMWYx8FhwmHkxMNkYsxT3kt9twOMh1iUi0iDKK6q47/lvWLx2N83Skrj7hlM5sWNzr8tqVOrSgyoGPuDQz4vqwyTgEWAF7mdi03AHZIiIeKK03MeUZ+aycuNemjdrwr0TT6VTdobXZTU6oQbUYGCotXZpfezUWntKtZ/LgZuDXyIiniosqeSup+ewdut+WmemcO9NI2jfJv3oL5R6F2pAWUB9WxGJaXsLy5n01Gw25xbRrlUa9944grYtNajYK6EG1APAC8aYJ4H1gK/6k9baD+u7MBGRcNq1t5Q7n5rNzvwSOrZtxj0Th9MqM/XoL5QGE2pAvRr8/vBhngvpRl0RkUi1bVcRf3xqDvn7y+jWIZMpPx1OZvp3Js2RMAspoKy1WthQRGLS+m37mfz0HAqKKzmpc0smTxhG09Qkr8sSjr7cRp0YY2odbi4iEklWbtzDHf/7NQXFlQwwWdx9w3CFUwSprQd1qzHmd8Bfgc+stb7DNQouange7ppQpcAX9V6liEg9W7A6j/tfmE+lz8+Ivjn8+sqBJCXqYlEkqW25jR8bYy4EHgQ6GWNm4t6rlI87YWwb3FV2hwNbgHustW82eMUiIsfpq0XbePSVhfgDDmcOOYGbf3IKCfE158EWr9X6GZS19h3gneBKuefghlFb3JkkcoEFwAPW2lkNXKeISL34cPZGpr69FMeBC0efyHXn9SIuTuEUiUIdJDGTo6yUKyISyRzH4fXP1vCPj1YDMO7cXowd093jqqQ2IS9YKCISrQIBh6ffW8b7/91IXBz8z8X9+OHwzl6XJUehgBKRmOarCvDYawv5atF2EhPiue3KgYzod9Rl7CQCKKBEJGaVVVTx4IvzWWh3kdokgTuuG0q/7m28LktCpIASkZi0v6iCKc/OZd3W/WSmJ3PXhOFaLiPKhBxQxpgsoC+QhDvM/CDNxScikSR3TwmTps1hZ34JbVumcfcNw8nRjORRJ6SAMsaMx10B93C3WGsuPhGJGBu2FzD56TnsL6qga04md/10GC0yUrwuS45BqD2o3wBPA7+31hY1YD0iIsds8Zpd3P/CfMoqquh7YmvuuG4IaSmauihahRpQHYHHFU4iEqlmLtjKY68twh9wGHlKe265vD9Jibq4E81CDahPgDOAtQ1Yi4hInTmOw9sz1vHCBysBd3aIa8/tRbymLop6oQbUEuBRY8z5wBqgsvqT1trb67swEZGj8Qccpr2zlA9nbyIuDsaf35sLRnbzuiypJ6EG1ChgHpCKO0FsdU69ViQiEoLyyioe/scC5q3IJSkxnl9dMYDv9WvvdVlSj0Kdi+/0hi5ERCRU+4squOe5uazZsp/01CTuvH4oJ3dt5XVZUs/qch9UW9w1n07GXehwFfC0tXZDA9UmIvIdW/OKmPLMXPL2lpLVIpW7fjqcjm2beV2WNICQVucyxgzB/ezpQtz1oHbjLlK41BgzqOHKExH5f8vX53P7E7PI21vKiR2b8/AvRiqcYlioPahHgFeBm6y1Bz9zMsY8CTwE6BKgiDSomQu38fhri6jyBxh6cjtuu3IgKU00W1ssC/XsDgImVA+noCdwFy0MiTHmPOB+oAuwC/iztfYpY0wy8CQwFvADj1prHwj1fUUkdjmOw2ufWF75xALwo9O6Mv783loBtxEI6RIfsBPofJjtXYGQbt41xmQDbwK/tdY2A34CPGaMGQBMAQzQDRgMjDPGXBNibSISo3xVfh59dSGvfGKJj4MbftyHG37cR+HUSITag3oZmGaMuQWYG9w2HPhL8LmjstbuNMa0sdYWGWPigVZAFW7AjQOutdbuA/YZYx4GJgIvhX4oIhJLCooreODF+azYsIeU5AR+c/UghvRq53VZEkahBtR9QA7wOm6vKw7w4V7iuyPUnQXDKQ0oCO77T7gDLrKBldWargb6hPq+IhJbtuYVcfezc8ndU0rLjBQmjR9Ktw5aKqOxCfU+qErgp8aY23AvxZUB66y1Zcewz3KgKe7SHR8G3wugtFqbUiDtGN5bRKLc4jW7ePDF+ZSUV9GtQyZ/vH4orTJTvS5LPHDEgDLGnAN8aq31BX+uqaMxBqjbelDW2gDuVEnfGmOm4Q7AAHeWigPSgOJQ31NEYsOHszfy1DvLCAQchvfJ5leXD9BIvUastjP/PtAOd7Td+7W0C2k9KGPMKNzReQOrbW4C7ANycXtm24Pbe3LoJT8RiWF+f4Cn31vOB19vBGDsmO5cffZJmvC1kTtiQFlr4w/383FYDLQ3xvwKeBwYCozHvfk3F5hsjFkKpAO3BduISIwrLvPxp5fms3jNbhIT4vn5JacwZlBHr8uSCBDqTBJfGGO+8wmlMaaNMSak+6CstQXAOcBFwF5gGu69VV8Ck4DlwApgPvAWMDWkIxCRqLVtVxG3Pf4li9fsJjM9mftvGqFwkoNq+wxqNNAr+HAUMNEYU/Oep5Nw710KibV2IfC9w2wvB24OfolII7BgdR4PvfwtJeVVdM7O4M7rh9K2pcZGyf+r7TOoPbiX2uKCXzfjzvJwgIM7kOHXDVadiMQcx3F476v1PP/vFQQcGN4nm1svH0CqBkNIDbV9BrUMd6YIjDEzgIuCN9KKiByTCp+fJ99YzMwF2wC4/CzDZWcaDYaQw6rtEl+atfbAvUnnHth2uLbV2omIHFb+/jLue+Eb1m3dT5PkBG69bAAj+uV4XZZEsNr61EXGmGxr7S7cS3mHWzk3jhCHmYtI47Viwx4efGk++4sqyGqZxp3XDaFLTqbXZUmEqy2gxuCOtgMtpyEix8BxHD6cvYmn312GP+DQ98TW3H71IDLTm3hdmkSB2j6D+vJwPwMEl8foC6yx1hY2XHkiEq0qfX6mvr2UT7/ZAsAFI7tx3Xm9SEioj9sqpTEIadiMMeZE4Fngt8BSYDZuQBUYY8621s6t7fUi0rjs3lfGAy9+w9qt+0lOSuDnP+nH6IG6v0nqJtRxnU/gLouxCbga6IA7NdF1wKPAqQ1RnIhEnyVrd/Pnl7+lsKSSrJZp/GHcYM1ELsck1L72acCt1tpc4MfAB9batcDTwCkNVZyIRA/HcXh7xlomPTWbwpJKBpgs/nLLKIWTHLNQe1DlQJIxpinurBLXB7e3w13bSUQasdJyH4+9tog5y3YCcMn3e3DFD3pq5Vs5LqEG1Me4vaUi3LWa/m2MOQN3QtfpDVSbiESBzTsLuf+Fb9iRX0JaSiK3Xj6AYb2zvS5LYkCol/gmAt/i9qTOtdaWAIOBmcAtDVOaiES6GQu28uu/fsWO/BI6Z2fwl1tHKZyk3oS6om4x8EsAY0yGMaa5tfbBBq1MRCJWpc/PtHeX8fHczQCcPrAD/zO2HynJmk9P6k/I/zcZY24C/gDkBB/vAh5XUIk0LjvzS3jwpfls2F5AUmI8Ey/sw1lDOxEXp8+bpH6Feh/UbcAfgfuA/+JOcTQC+J0xpsxaq8UFRRqBr5fu4K//WkRpeRXtWqXxu2s0hFwaTqg9qJuBG621r1bb9rUxZjNwL1r9ViSm+ar8PP/+Sv49awPgLpHxi0v7k56a5HFlEstCDag2uCvd1rQA96ZdEYlRuXtK+PPL37J2634SE+K47ryT+dFpXXVJTxpcqAG1HPgJ8ECN7ZcCq+u1IhGJGLMWb+fJNxZTWl5FVotUfnvNYHqc0MLrsqSRCDWgJgEfGGOGA3OC24YDPwQuaojCRMQ7FT4/z763nP/M2QQEL+ldcgrpacme1iWNS6jDzD8J3pj7c9y5+MqAVcBga+2SBqxPRMJsc24hD738LZtzi0hMiGfCBb0559TOuqQnYRfyMHNr7VfAVw1Yi4h4yHEcPp67maffW06lz0/7Nk25/erBdG2vhQXFG7Uu+Q48BowFKoB3gN9p/SeR2FNUWskTry8+OJfe9wefwA0X9iG1iW68Fe/U9n/fFOBHwJ8BP/AzoBXuwAgRiRHL1uXzyCsL2FNQTmqTRG4e249RAzQ4V7xXW0CNBa6w1s4AMMZ8CXxljEmy1vqOZWfGmDOBB4HuwC7gIWvtU8EVep8M7tMPPGqtrTliUETqka8qwKufrObNL9biONCzUwt+feVA2rVq6nVpIkDtAdWBQ4eQz8edXLYtsK2uOzLGdATeAsYB7wEDgY+NMZuA0bgLIHYDMoGPjDHbrbUv1XU/InJ0W/OKePSVBazbVkB8HFx6puGyM3toOXaJKLUFVAJubwYAa61jjKkAjnWcaWfgFWvtO8HH840xM3GnTBoHXGut3QfsM8Y8jDuDugJKpB45jsN/5mzi2ekrqPT5yWqZxq8uH8DJXVt5XZrId4TtE1Br7Sxg1oHHxpiWuCv1vgxkAyurNV8N9AlXbSKNwd7Ccv76r0UsWL0LgDGDOnLDj/vQVNMVSYQ6WkBda4wprtH+KmNMfvVG1tq/12WnxphM3IUO5+FOlwTuQohU+zmtLu8pIkf29ZId/O3NJRSVVpKemsT/XNyP0/q397oskVrVFlBbgJtqbMsFrquxzQFCDihjTA/cz6BWAlcCqcGnUqs1SwOKEZHjUlxayVPvLmPmAvdj4/492vDLy/rTKjP1KK8U8d4RA8pa27m+d2aMGYkbTlOBP1hrHaDcGJOLO0hie7BpTw695CcidbRw9S7++voi9hSUk5yUwPXn9eKcEV00I4REjbB9BmWM6Qa8D9xhrX2ixtMvA5ONMUuBdOA2tISHyDEpLffx/Psr+WjOJgBMpxbcevkA2rdJ97QukboK523iNwPNgAeMMdXvcfob7mS0jwArcIeyT8PtZYlIHSxZu5u//msRu/aVkZgQxxU/6MlFo0/U8HGJSuEcxfcr4Fe1NLk5+CUidVRWUcWLH6zkg683AtC1fSa3XNafLjmaR0+ilybaEolyi9fs4ok3lrBrbymJCXFceqZh7JjuJKrXJFFOASUSpUrKfDz//go+nrsZUK9JYo8CSiQKzVu+k/99eyl7CspJTIjn8rMMF51+onpNElMUUCJRZF9ROdPeWcZ/l+wAwJzQgp9fegqd2mV4XJlI/VNAiUQBx3H4fP4Wnp2+guIyH02SE7j67JM473tdSYjXfU0SmxRQIhFux+5i/vbmEpauc2cYO6VHG24e20/LYkjMU0CJRChflZ+3Z6zjX5+twVcVIKNpMhMu6M3oAR00G4Q0CgookQi0bH0+f39zCdt2uVNSjhnUket/dDKZ6U08rkwkfBRQIhGkoLiC599fwefztwLQvk1Tbrq4H/26t/G4MpHwU0CJRIBAwOHjeZt56YOVFJf5SEyI55IzujP2jO4kJSZ4XZ6IJxRQIh5bu3UfU99eypot+wF3SYwbL+pLjiZ3lUZOASXikcKSSl76cCWfzNuM40CrzBQmXNCbEX1zNAhCBAWUSNj5Aw4fz93EP/6ziqJSHwnxcVwwqhuXntmDtBQtvy5ygAJKJIyWr89n2rvL2LijEIB+3Vsz8cK+dGzbzOPKRCKPAkokDHbtLeWFD1Yya7G7aHSbFqmM/1FvTu2brct5IkeggBJpQOUVVbz5xVrembmOyqoAyYnxjB3TnQtPP5GUZP3zE6mN/oWINIBAwGHGgq289OEq9haWAzDylPaMO68XWS3SPK5OJDoooETq2dJ1u3l2+go2bC8A4MSOzfnpBb3p1aWVx5WJRBcFlEg92ZJbyIsfrOKblbkAtM5M4ZpzezGqfwfiNeO4SJ0poESO097Ccl75eDWfzttMwIGU5ATGjunOBaO66XMmkeOgfz0ix6ikzMdbM9YyfdYGKir9xMfHcfbwTlx+lqFFsxSvyxOJegookTqq9Pn5cPZGXv9sDUWlPgCG9W7HNef00v1MIvVIASUSoip/gM/nb+G1Tyz5Be7IvJO7tuLac3vRs3NLj6sTiT2eBJQxZgjwvrU2K/g4GXgSGAv4gUettQ94UZtITf6Aw6zF23n149XsyC8BoHN2BteccxKDTmqrG21FGkhYA8oYEweMBx6u8dQUwADdgLBUdPcAAA9zSURBVEzgI2PMdmvtS+GsT6S6QMBhzrKd/PPj1WzNKwIgu3VTrvphT77Xr71G5ok0sHD3oKYA5wL3AndW2z4OuNZauw/YZ4x5GJgIKKAk7AIBh3krdvLqJ/bgnHlZLVK55PuGMwZ3JDEh3uMKRRqHcAfUVGvtJGPM6AMbjDHNgWxgZbV2q4E+Ya5NGrnDBVOrzBQu+X4PzhzSiaREBZNIOIU1oKy1Ow6z+cCqbKXVtpUCmg9GwsIfcJi9ZAf/+syyOde9lNcqM4WfjOnOmUM7kZykFW1FvBAJo/hKgt9Tq21LA4o9qEUakSp/gJkLtvHmF2vZvtv93611ZgoXj+nOWQomEc95HlDW2n3GmFzcQRLbg5t7cuglP5F6U+Hz8+m8zbw9cx2795UBkNUyjZ+M6c4ZgzuSlKhgEokEngdU0MvAZGPMUtxLfrcBj3tbksSaotJKPvx6I//+7wYKiisB6JCVztgx3Rk1oIMGP4hEmEgJqEnAI8AKIB6YBkz1tCKJGXl7S5n+1Xo+mbeZ8ko/4M4wfskZ3Rl6craGi4tEKE8Cylo7E2he7XE5cHPwS6RerNmyj3dmrmP20h0EHHfbgJ5ZXHz6ifTp1lo32IpEuEjpQYnUC78/wNzlubz31XpWbdoLQEJ8HKcPaM+Fo0+kS06mxxWKSKgUUBITikor+XTeZt7/euPBgQ9NUxL54fDOnPe9rrRunnqUdxCRSKOAkqi2cUcB/561gS8XbqOyKgBATuumnH9aV8YMPoHUJvpfXCRa6V+vRB1flZ+vl+7kw683HryMB9C/RxvOO60rg3q21cAHkRiggJKosSO/mE/mbuaz+VsODhNPbZLIGYM7cu6ILnTI0lpMIrFEASURzVflZ+7yXD6eu4kla/MPbu+Sk8E5p3Zh1IAOuownEqP0L1si0qadhXw6bzMzFmyjqNTtLSUnJXDaKTn8cFhnTKcWGiYuEuMUUBIxCksq+WrRNj6fv4V12woObu+ak8mZQ09g9MCOpKcmeVihiISTAko85avyM39lHjMXbmP+ylyq/O4dtU1TEhk5oANnDe3EiR2aH+VdRCQWKaAk7AIBhxUb9/Dlwm18vWQHxWU+AOLj3Jkevj/oBIb2bqfZxEUaOQWUhIXjOKzdup9Zi7cza/F29hSUH3yuS04Gpw/syMj+7WmVqRtqRcSlgJIG4zgO67bt5+slO/jvkh3k7f3/NSmzWqQysn8HRg3oQOfsDA+rFJFIpYCSeuUPOKzetJc5y3YyZ9kOdgWnHQJomdGEU/vmMKp/B43CE5GjUkDJcavw+VmyZjdzl+/km5W5B2+iBTeUhvfJ4Xv9cujVpZVmeBCRkCmg5Jjk7y/j21V5fLMylyVr86n0+Q8+165VGsN6Z3NqnxxMpxYKJRE5JgooCYmvKsDqTXtZsDqPBat3sWln4SHPn9ghk6G9sxnWO5tO7Zrp8p2IHDcFlByW4zhszSti8drdLLK7Wb4+/+BqtAApyQn0696Gwb3aMeikLI2+E5F6p4ASwA2knfklLFu/h2Xr8lm6bjf7iioOadOpXTP6mywG9WxLr64tSUrUfUoi0nAUUI1UIOCwObeQlRv3snLDHpZvyGdv4aGB1KJZE/p1b8MpPdwv9ZJEJJwUUI1EcZmPNVv2YTftZfXmfdjNeykprzqkTUbTZHp3a0Wfbq3p170NHbLS9VmSiHhGARWDyiqq2LC9gPXb97N2637WbtnH9t0l32mX1SKVXl1a0atLS3p1bcUJbTW4QUQihwIqijmOw97CcjbuKGTjjgI27Shkw44Ctu8uxnEObZuUGE/X9pn07NSSnp1b0LNTS1o31yU7EYlcCqgoEAg47CkoZ+uuIrbtKmJbXjFb8orYvLPw4ESr1SXEx9EpO4NuHTLp3rE53Tu2oFN2BkmJ8R5ULyJybCImoIwx/YCpQF9gA3C9tXa+t1WFT5U/QP7+MvL2lJK7t4TcPaXsyC9mx+4SduSXHHIjbHXpqUl0ys6gS04GXXIy6ZKTQefsDI2wE5GoFxEBZYxJBt4DHgNGAhcDnxhjOllrC2t9cYRzHIeyiir2F1Wwt7CcfcHv+fvL2FPgft+9r5S9heUEnCO/T2Z6Mh2ymtEhK50OWc3o1K4ZJ7RrRsuMFH1uJCIxKSICChgNJFlrHws+fs0Y8zPgUuDpcBTg9wcIOG6gBBwHv9/BH3Dw+wP4qgJUVvnxVQWo8Pkpr6iivNL9XlpRRUmZj9Jy93tRaSXFpT4KSyopKKmgsKQSX1XgqPuPj4NWmSm0a9WUti3TyG7dlHYt08hpk05Om3StJCsijU6kBFQvYFWNbauBPuHY+XP/XsE7M9c12Ps3SU6gRbMmtGiWQouMJrRslkKr5qm0zkyhVWYqbVqk0rp5KokJ+oxIROSASAmodKC0xrZSIC0cO09MiCMhPo64OIiLiyMOSEiIJyE+jsSEeBIT42mSFE9SYgJNkhJokpxASnIiKckJpKUm0TQlkdSURNJTk2mWlkR6ajLpaUk0T29CRnoyKcmR8p9ZRCR6RMpvzhKg5pjnNKA4HDu/5pxeXHNOr3DsSkREQhQp15RWAqbGtp7B7SIi0ghFSg9qBhBnjLkVeBJ3FF9f4B1PqxIREc9ERA/KWlsJnI0bTHuBO4AfW2t3e1qYiIh4JlJ6UFhrlwPf87oOERGJDBHRgxIREalJASUiIhFJASUiIhEpYj6DOg4JALm5uV7XISIidVTtd/d3ZriOhYDKBrjyyiu9rkNERI5dNrC++oZYCKj5wGnATuDwa1KIiEikSsANp+8srxTn1Fx6VUREJAJokISIiEQkBZSIiEQkBZSIiEQkBZSIiEQkBZSIiEQkBZSIiEQkBZSIiEQkBZSIiESkWJhJIiTGmH7AVNyVejcA11trv3PnsjHmBOBZYBiwC/i5tfbDcNZ6LOpwfGOAT4Gyapv/ZK29JyyFHidjzBDgfWtt1hGej8rzd0AIxxd1588YcybwINAd95w8ZK196jDtovLc1eH4ou7cARhjzgPuB7rgHt+fw3X+GkUPyhiTDLwH/AtoDtwHfGKMyThM89eApUAr4KfAa8aYruGq9VjU8fgGAG9Ya9OrfUX0PxAAY0ycMWYC8AmQXEvTqDt/UKfji6rzZ4zpCLwF3Iv7/+blwAPGmB8cpnnUnbs6Hl9UnTsAY0w28CbwW2ttM+AnwGPGmAGHaV7v569RBBQwGkiy1j5mrfVZa18DVgCXVm9kjOkBDAImWWsrrbVfANOB8eEuuI5GE8LxBQ0EFoezuHoyBbgJ9xfBYUXx+YMQji8o2s5fZ+AVa+071tpAsFc/ExhRvVEUn7vOhHB8QdF27rDW7gTaWGv/Y4yJxw2fKqCoeruGOn+N5RJfL2BVjW2rgT6HabfFWltSo92QBqytPoR6fOD+FdfGGHMTEIfb67rTWlvRsCUet6nW2knGmNG1tInW8wehHR9E2fmz1s4CZh14bIxpiTu588s1mkbluavD8UGUnbsDrLVFxpg0oAA3M/5krV1bo1mDnL/G0oNKB0prbCsF0o6xXaQJqW5jTCKwDXgHOAkYA3wfiOjLDADW2h0hNIvW8xfS8UXz+QMwxmTi/lU9D/eSdHVRe+4OqO34ov3cAeVAU2AwcL0xpmbPqEHOX2PpQZUAqTW2pQHFx9gu0oRUt7W2Cjij2qZ1xpj7gD8BtzdoheERrecvJNF8/oKXgN4DVgJXWmsDNZpE9bk72vFF87kDCB5PJfCtMWYacAHugIgDGuT8NZYe1ErA1NjWM7i9ZrsTjDGpR2kXaUI6PmNMe2PMw8FBFQck4/51FAui9fyFJFrPnzFmJG6v4l1grLX2cPVG7bkL5fii+NyNMsYsqLG5CbC/xrYGOX+NpQc1A4gzxtwKPAlcjDsc+53qjay11hizBLjPGPN74FTcvxSGh7neugrp+IA9wJVAqTHmbtxho3cCz4Wx1gYTxecvVFF3/owx3YD3gTustU8cqV20nrtQj48oPHdBi4H2xphfAY8DQ3EHPlxYvVFDnb9G0YOy1lYCZ+P+4t4L3AH82Fq72xhzpTGmejf0YtxrxLuAZ4Dx1trl4a65LkI9vuBfdmcDI3H/wXwFvAE86knh9SAWzl9tYuD83Qw0wx16XVzt608xcu5COr4oPXdYawuAc4CLcH+3TAMmWGu/DMf504q6IiISkRpFD0pERKKPAkpERCKSAkpERCKSAkpERCKSAkpERCKSAkpERCJSY7lRV6TeGWNeAMbV0mQK7szWM4Bm1tqwTNtjjEkAvgausdauqaVdPDAXuNpaa8NRm0hdqAclcux+CWQHv0YHtw2ptu1hYHbw55LDvL6h/AJYUls4wcH51e7GXehSJOLoRl2RemCM6Q0sA7pYazd5WEcKsAUYE+pd/MaY9bh3/c9syNpE6kqX+EQaUHB9p4OX+IwxDu6qq7/HneD3W+Aq4DfA1UAh8Htr7cvB1zcDHgHGAg7wBfDLWpbnuAzYXz2cjDF/BG4A2uCuG/YHa+1/qr3mHdze4Mx6OGSReqNLfCLh9yBwCzAMOAFYiBtMg4G3gaeMMenBttNwg+wHwCjckPo4uL7Q4ZwLfHTggTHmwuC+rsKdXfoD4A1jTEa113wEfL+W9xTxhAJKJPz+Zq2dYa1djDsTdjFur8biTh6aCnQxxnTF7RFdYa2dH+wVXY27zPgPj/Deg4AV1R53BiqAzcFLj3fjTvzpq9ZmJe6Ccz3r5ehE6on+YhIJv3XVfi4FNllrD3wYfGB9oCZAp+DP1phDlvtKw+1VvX+Y924L5Fd7/A/ckYYbguv6TAeet9aWVWuzJ/g9q47HIdKg1IMSCT9fjcc1V5c9IDHYtj9wSrWvHsDzR3hNAIg78MBauxsYiNvjmg1cCywNDuo44MDvAX/IRyASBgookci1CkgCmlpr11lr1wE7gYdwQ+pwcnEHQwBgjLkImGit/cRa+0vcnlcR7ho/B7Sp9lqRiKFLfCIRKrhK6XTgJWPMzcBu4D7cwRWrj/CyBUC/ao8TgIeMMXm4IwaHAe2CPx/QD9jHoZceRTynHpRIZBuHGybvAvOBTOBMa+3+I7T/AHe0HwDW2jeAybi9rjXAvcDPrLVfVHvNSOAja60u8UlE0Y26IjHEGJMGbAJ+aK1dGEL7eGAz7kjBWQ1cnkidqAclEkOstaW4vaWbQ3zJBcAGhZNEIgWUSOz5C9DX1BibXlOw93QHcGNYqhKpI13iExGRiKQelIiIRCQFlIiIRCQFlIiIRCQFlIiIRCQFlIiIRKT/A5LSCQDPSVg8AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "def plot_position(results):\n",
    "    plot(results.x, label='x')\n",
    "    decorate(xlabel='Time (s)',\n",
    "             ylabel='Position (m)')\n",
    "    \n",
    "plot_position(results)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plotting `v`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO3dd3xUdb7/8VcKJaFKr9L5hhKKBRUVsaCAtAiKiq7Y137v7t772133rte7xS2uq2tZC1bsIEVAEbvYQaWGfOjd0HsSUmZ+f5xBxwhxApk5M5P38/HIg8zJyczneGLe+Z75nu8nJRgMIiIiEm9S/S5ARETkcBRQIiISlxRQIiISlxRQIiISl9L9LuBYOedqAScD3wJlPpcjIiKVkwa0BOaZ2cHwLyR8QOGF01y/ixARkWNyJvBx+IZkCKhvAV544QVatGjhdy0iIlIJ+fn5jBs3DkK/y8MlQ0CVAbRo0YI2bdr4XYuIiBydH71Fo0kSIiISlxRQIiISlxRQIiISlxRQIiISl2I+ScI51xBYBPzezJ5xztUEHgLG4L1Jdp+Z3VOVr7l37162bt1KSUlJVT6tb+rUqUObNm1ITdXfFyKSvPyYxfco0Drs8d2AAzoBDYDZzrlNZvZcVbzY3r172bJlC61btyYjI4OUlJSqeFrfBAIBNm3axPbt22nWrJnf5YiIRE1M/wR3zl0F1AcWh22+CviTme0ys7XAvcCNVfWaW7dupXXr1mRmZiZ8OAGkpqbSvHlz9uzZ43cpIlLNbdq2n70HiqP2/DEbQTnnOgB3Af2B2aFtDfGWuMgN2zUPyK6q1y0pKSEjI6Oqni4u1KhRg9LSUr/LEJFqatO2/Tw7K5fPFn9Ldqcm/Pnm06PyOjEJKOdcGvA88Cszy3fOHfpS3dC/BWG7FwCZVfn6yTByCpdsxyMiiWHP/oO8NMeY/dlaygJBatVM45yTordAQqxGUP8DmJlNKbf9QOjf8CFOJrA/JlWJiMhPOlhSxusfrWLSuysoPFhKagqcf0o7Lr/A0bhB9K5QxSqgLgVaOecuCj2uBzwC9APy8SZJbAp9LYsfXvITEREfBAJBPvh6AxPfWMb2PUUAnJjVjPHDetC+Zf2ov35MAsrMssIfO+cWAPeHppnvB+5yzi3Cu+T3K+CBWNQlIiKHt3D5Np6asZTVm70JWR1bNeCa4T3o3bVpzGqIh8Vifw/8A1iKN6vwcbyp6NXKL3/5Sxo1asSdd94JQFlZGWeeeSb//Oc/OeWUU3yuTkSqi/X5e3l6Zi7zl20BoEmD2lwxpBtnn9iW1NTYvv/tS0CZWZ+wz4uAW0IfUXf3hM+/+w8fCyd1a85d1536k/vl5OTw61//ml//+tekpaXxySefULt2bfr16xeDKkWkutu1r4gX3zLmfL6WQBAyaqVz8bldGDGgE7VqpPlSUzyMoATo378/qampfPHFF/Tv358ZM2YwfPhwzdgTkagqKi5l+oereO39FRQeLCM1NYWhp7XjsvOzaFivlq+1VbuAimQ044fU1FRGjBjBjBkz6Nu3L++88w5TppSf9CgiUjUCgSDvf7WBiW8uY0doAsQpPVpw1YXdadu8ns/VeapdQMWznJwcLr30Uk499VS6du1Khw4d/C5JRJLQwhWhCRCbvAkQndp4EyB6dY7dBIhIKKDiSKdOnWjXrh33338/119/vd/liEiS2bBlH0/PXMq83O8nQFw5tDsDT2gT8wkQkVBAxZmcnBzuuecehg4d6ncpIpIk9uw/yItv5TH783UEAkEyaqUx5pyujDzLvwkQkVBAxZlx48Yxbtw4v8sQkSRQXFLG63NXM+nd5RQUeStAXHBqO8YNzuK4erX9Lu8nKaBERJJMMBhk7oJNPDsrl627CgFvBYirh/egXYvorwBRVRRQIiJJJG/tTia8vgRbtwuAdi3qcc2InpzgEq9/nAJKRCQJbNlZwLOzcpm7wFvWtGG9WlwxuBvn9TuetDicABEJBZSISAI7UFjCpHeX8/rc1ZSUBqiZnsqogZ0ZfXZnMmvX8Lu8Y1ItAioYDCbVigzBYNDvEkTEZ2VlAeZ8sY4X3spjz36vq+3AE9vwsyHdaXpccjRpTfqAqlGjBoWFhWRmVmkPRF+VlJSQnp70p05EjuCrvC08+fpSNmzZB0D3Do24dkRPuh5/nM+VVa2k/y3XrFkzNm3aROvWrcnIyEj4kVQgEGDLli00aNDA71JEJMbW5e/lqdeX8rVtBaBF40zGD+tB/+yWCf+77XCSPqDq1/emVG7evJmSkhKfq6kaderUoUmTJn6XISIxsmf/QV6YncdboZXGM2unM/a8rgw/syM10uP3RttjlfQBBV5IHQoqEZFEUVJaxoy5q3nlndCNtqGVxi+/IIsGdf1daTwWqkVAiYgkkmAwyKeLv+WZmUvJ31EAwAlZzbgmwW60PVYKKBGROLJy424mTF/C0tU7AGjbvB7XjujBiVnNfa4s9hRQIiJxYOfeIia+sYx3568nGIR6mTUZNziLwae2Iy0t1e/yfKGAEhHx0cGSMqZ9uJLJ766gqLiM9LQUhp3RkbGDHHUzEvtG22OlgBIR8UEwGOTjBZt5etZStoUWdD21ZwuuHtaDVk3r+lxdfFBAiYjE2PL1u5gwfQnL1u4EoEOr+lw3smfcdbT1mwJKRCRGduwp5Lk3lvHe/A0ANKxbiyuGJPaCrtGkgBIRibKi4lKmfbiKye+t4GBxGelpqYwc0JFLzuua8Au6RpMCSkQkSg41Dnx6Zi7bd3vvM/Xv1ZKrh/WgReM6PlcX/xRQIiJRcLj3ma4fmU12Zy1TFikFlIhIFdq5t4hnZ+XqfaYqoIASEakCxSVlTP9oFZPeXU7hQe9+phFnduKS87pSp5rfz3S0FFAiIsfg0Lp5T89Yypad3rp5p/ZswdXDe9Cqie5nOhYKKBGRo7Rm8x6emLaExau2A9C+ZX2uG9GT3l11P1NVUECJiFTSnv0HeX52HnNC/ZnqZdbkiiFZXHBK9V03LxoUUCIiESopDTDrkzW8PCePA6H+TMPP6MBl5zvqZdb0u7yko4ASEYnA/GVbmDB9CZu27QfgBNeMa0f04Phq1J8p1hRQIiIV2LRtPxOmL2H+si0AtGpSh2tH9uTkbs1JSdG08WhSQImIHMaBwhJeftuYMXc1ZYEgmbXTuXSQY9gZHamRrveZYkEBJSISJhAI8s689Ux8Yxm79x8kJQXOP6UdVwzJ4rh6tf0ur1pRQImIhOSu2cET0xazcuMeALq1b8QNo7Lp3Lahz5VVTwooEan2duwp5OkZuXz4zUYAGjeozfhhPTirb2u9z+QjBZSIVFvFJWVM+9BbnqiouIwa6alcdHZnxpzdhdq19OvRbzoDIlLtBINBPl+Sz1MzlpC/w1ue6LTsllwzXG0w4okCSkSqlfX5e3li+hIWLN8GwPEt6nHDqGx6d9HyRPFGASUi1cL+whJempPHzI/XEAgEqZNRgysGZzHktPZanihOKaBEJKkdmjb+3Bu57NlfTGoKDDmtPeMGZ9Ggbi2/y5MKxDSgnHPDgD8DHYCtwN/M7DHnXE3gIWAMUAbcZ2b3xLI2EUk+eWt38ti0xazcsBuAHh0bc8OobDq2buBzZRKJmAWUc64lMBnIMbM3nXMnAJ845+YBFwMO6AQ0AGY75zaZ2XOxqk9Ekkf5rraNG9Tm6mE9GKBp4wklZgFlZt8655qa2T7nXCrQGCgF9gFXAePNbBewyzl3L3AjoIASkYiVlAaYMXc1L79tFB4sJT0tlZyBnbj43K5kaNp4wonpGQuFUyawJ/TafwW2AS2B3LBd84DsWNYmIont67ytPD5t8Xerjffr3oLrRvakZRNNG09UfvxJUQTUAXoBbwCFoe0FYfsUAJkxrktEElD+jgNMmL6EL5bmA9C6aR2uH5XNiVnNfa5MjlXMA8rMAkAxMN859zhwUuhLGWG7ZQL7Y12biCSOouJSJr+7gikfrKSkNEBGrTQuHeQYfmYnrTaeJGI5SeIsvNl5J4ZtrgXsAvLxJklsCm3P4oeX/EREAG8ViE8WbebJ15eyfbd3AebsE9swflgPGtXXauPJJJYjqAVAa+fcL4AHgFOAa4EcvIC6yzm3CKgL/Cq0j4jId9bl7+XxqYtZtHI7AB1bN+DGnGy6d2jsc2USDbGcxbfHOTcU+BdwF7ABuM7MPnTOfQH8A1gKpAKPA4/GqjYRiW8HCkt4aY4x4+PVBAJB6mXW4Mqh3Tn/lHakpWraeLKK9Sy+r4EzDrO9CLgl9CEiAnirQLz/1QaemZXL7n1e88Ah/dtzxeBu1K9T0+/yJMp0Y4CIxKWVG3fz2JRF5K3bBXjNA2/MyaZTGzUPrC4UUCISV/YeKOb5N5cx+/O1BIPQsF4trh7Wg7NPbKNVIKoZBZSIxIWyQJC3v1jHc28sY19BMWmpKQwf0JHLzndk1q7hd3nig4gCyjmXDQzBu2epGd6CrvnAPGCmma2MWoUikvTy1u3ksSmLWLlxDwC9Ojfhxpxsjm9R3+fKxE8VBpRzbgBwN9Af+BLv3qSVQBrQBLgC+Jtz7kPgD2b2UXTLFZFksmf/QZ6dlcvbX64HoEmD2lw7sien92qly3ly5IByzj0F9MBrg5FjZruPsF994DLgfufcIjMbH41CRSR5lJUFePOztTw/O48DhSWkp6WQM7Azl5zbldpa1FVCKvpJmGVm1/zUE5jZXuAx4DHn3Jgqq0xEklLumh08NmUxqzd7l/NOcM24ISeb1k3r+lyZxJsjBpSZvVbZJzOzycdWjogkq137inhm5vc9mpodl8F1I7M5tWcLXc6Tw4p0kkQm8P+A581shXPuMWAc3vtS48zs2yjWKCIJrKwswKxP1/DC7DwKikqpkZ7KRWd3Zsw5XahdU5fz5Mgi/el4ADgHeNU5NxKvweBtwCi8pYsujk55IpLIlq7ewaNTFrH2270AnNStOTeMylaPJolIpAE1EhhmZkudc78B3jazJ5xznwCfRa88EUlEP7qc1yiTG0dl069HC58rk0QSaUBlAFtCrdovAP43tD2Id0+UiMhhL+eNPrsLY87tQq0aaX6XJwkm0oCah/ce1DbgOGCqc64V8Afg8yjVJiIJJHeNdzlvzWZdzpOqEWlA3Qq8CLQHbjazzc65B/GaDOZEqTYRSQC79x3kmVlLeXfe97PzbghdztPsPDkWFd2o2x/43MwCZpYL9Cm3y2/NbF9UqxORuFUWCDL7s7VMfHNZ6GbbVEafo9l5UnUq+im6F3DOubnAHGBO+Jp7CieR6svW7eTfUxaxauP3N9vemJNNK91sK1Wooht1+4eWMToXGAT8wjmXBryNF1jvmtmu2JQpIvFg74FinnsjlzlfrCMYhCYNM7h+ZE9Oy26py3lS5Soch4eWMZoa+sA51wE4H7gUeNQ5twpvyvnvol2oiPgnEAjy9pfreXZW7netMEYN7MSlg5zWzpOoqdRPlpmt4ft191KBk/FGVyKSpFZu3M2jry3C1nsXTHp1bsLPL+pF2+b1fK5Mkl3EAeWcG4i3unmtcl8qqMqCRCQ+7C8s4YU3l/HGp2sIBKFR/VpcO6InZ/Zprct5EhORrsV3P97SRuuBonJfDgL3VXFdIuKTYDDIh19v5MkZS9m97yCpqSmMPLMjl1+gzrYSW5GOoH4GXGNmz0azGBHx1/r8vTw6ZTGLV20HoFv7Rtw0uhcdWjXwuTKpjiINqAK8lctFJAkVHSzl5beNaR+uoiwQpH6dmlw9rAfnnNSW1FRdzhN/RBpQfwTudc7dGpooISJJIBgM8vmSfJ6YvphtuwpJSYEhp7XnyqHdqJdZ0+/ypJqLNKCWAX8GVjrnfvRFM9MqkCIJJn/HAR6bupj5y7YA0KlNA24e3Zuuxx/nc2UinkgD6nG8RWGfRrP2RBJaSWkZUz5YyatvL6e4NEBm7XSuHNKNIf07kKbLeRJHIg2otsAQM1sdzWJEJLoWrtjGv19bxKZt+wE4q28brh3Rg+Pq1/a5MpEfizSg3gYGAAookQS0a28RT76+lA+/2QhA66Z1uWl0L3p3aepzZSJHFmlAfQE85JwbDawESsK/aGb/XdWFicixKwsEmf3pGm/F8aJSaqanMnaQI2dgJ2qk661jiW+RBtQgvKaFdflx241glVYkIlVixYZdPPLaIlZu2A14DQRvzMmmRWM1EJTEEFFAmdnZ0S5ERKrGgcISJoaWKAoGoUmD2lw/KlsrjkvCqahh4V3A38ysMJIncs7VA/7LzH5fVcWJSOSCwSBzF2xiwvQl7Dq0RNGAjlx+QRYZWnFcElBFP7V7gKXOucnAFDP7vPwOzrkU4CTgCuAi4J9RqVJEKrR5237+PWURC5ZvAyCr3XHcPKa3liiShFZRw8L7Q+H038Ac51wp3g2724EUoCne6uYpwDPA6Wa2PuoVi8h3SkrLmPzuCia9t4KS0gB1M2owflh3BvVrpyWKJOH9VMPCjcDtzrnfAAOBE4HmQABvZt/dwPtmdjDKdYpIOQuXb+PfUxayadsBAM45qS3XDO9Bg7rlO+KIJKZIJ0kcAGaFPkTER7v2FfHk9O/vaWrTrC43j+5NducmPlcmUrX0zqlIgggEgrz1+VqenZX73T1NlwzqykUDu1AjPdXv8kSqnAJKJAGs2byHhycvxNZ5bddPyGrGTRf10j1NktQUUCJxrOhgKS/OMaZ/tIpAIEij+rW4bmQ2Z/RupXuaJOlF2vK9vZmtjXItIhLmy6X5PDp10Xd9moad3oErhnSjTobarkv1EOkIaqVz7nPgeWCSme2IYk0i1dr23YU8Pm0xny3+FoCOrRtwyxj1aZLqJ9KA6ghcDtwEPOCcmwO8AEyPdKUJAOfcIOAvQBdgK/B3M3vMOVcTeAgYA5QB95nZPZEfhkjiKwsEmfXxap6fvYzCg2Vk1Epj3OBuDDu9A2lpmgQh1U+k08zX4wXLX5xzPYFLgd8AjzvnpgITzeydip7DOdcWeA24CpiOd0/VW865tXj3WDmgE9AAmO2c22Rmzx3FMYkknJUbdvPw5AWs3LgHgFN7tuCGUb1oelyGz5WJ+OdoJklsBFbh9Ybqgre6+SDn3D5gvJl9doTvaw+8aGZTQ4/nOec+AE7HC63xZrYL2OWcuxe4EVBASVIrKCrhhdl5zPx4NYEgNGmYwc9zsjmlZ0u/SxPxXaSTJOoAI/FGTucDW4AXgd+Z2VLnXCrwCPAKcPzhnsPM5gJzw56zEXAmMBFoCeSG7Z4HZFf2YEQSyWeLv+WxqYvYsaeI1NQURmlhV5EfiPT/hK1AMTAFGGxmH4R/0cwCofelzozkyZxzDYDX8ZZL+iq0uSBslwIgM8LaRBLKtl2FPDZ1EV8szQegS9uG3DKmN53aNPS5MpH4EmlAjQdeP9yae865Zma21cym4AVYhZxzXfHeg8oFxgGHLrKHX2zPBPZHWJtIQvjxJIh0rhzSjaGndyBNC7uK/EikAfUy0ALYFr7ROXc8XtDUjeRJnHMD8MLpUeC3ZhYEipxz+XiTJDaFds3ih5f8RBLayo27eXjS95MgTstuyQ2jsmnSUJMgRI6kooaFlwE5oYcpwATnXPkRVDtgZyQv5JzrBMwE7jSzB8t9eSJwl3NuEV7Y/Qp4IJLnFYlnhQdLeWF2HjPmrtIkCJFKqmgE9TYwCC+cAApDH4cE8d5DeibC17oFqAfc45wLv8fpYeD3wD+ApUAq8DjeKEskYX2Zm8+jU7yVIFJTYMSAjlwxuJsmQYhEqKKGhduBawBC9yr93cwKjrT/TzGzXwC/qGCXW0IfIglt594iHp+2mE8Wbga8lSBuu7gPndtqEoRIZVR0iW8o8LaZlQDzgIHOucPua2ZvRKc8kcRRvh1G7ZreShDDz9BKECJHo6JrDTPxJkZsDX1+JEEgrSqLEkk06/L38vCkhSxb670le3L35vw8pxfNGuluCZGjVdElvtTDfS4i3ysuKePVd5bz2vsrKC0Lcly9WtyQk83pvdQOQ+RYRfxurXPuGmCvmU0OPX4VmGFmE6NVnEg8W7RyGw9PWsjm7QcAGHJae352YXfqqh2GSJWIdKmjO4Ff4q1mfshi4H7nXGMzuz8axYnEo70Hinl6xlLembcegLbN63LLmD706NjY58pEkkukI6gbgUvNbM6hDWb2B+fcfLw1+BRQkvSCwSAffrOJCdMXs2d/MelpqYwd1JXRZ3emRrrehhWpapEG1HHAusNsXwU0r7pyROLTlp0FPPLaQr7O2wpAj46NufXi3rRpVs/nykSSV6QB9Tnwa+fc9WZWCuCcS8O77DcvWsWJ+K2sLMDrc1fzwlt5HCwuo05GDa4Z3oPzTj6eVK2fJxJVkQbUr4B3gfWh5YiCeO0w0oEhUapNxFerNu7mobD1887o3YobRmVzXP3aPlcmUj1E2lF3ofPu0r0U6IbXemM68IKZ7YtifSIxV1RcyktvGdM+WkUgEKRJwwxuHt2Lk7u38Ls0kWol4mnmZrbDOfcWsAFvvbw8hZMkmwXLt/Lw5IXk7yggJQVGnNmRK4Zo/TwRP0Q6zbwu8CQwBijBW0A23Tn3NjDazA5Er0SR6Nt7oJgnX1/Ce/M3ANC+ZX1uu6QPXY8/zufKRKqvSP8svA/vPafT+H5SxMnABOCvwK1VX5pI9JWfOl4jPZXLznfkDOxMutbPE/FVpAF1EZBjZl+GbfvSOXcLMBkFlCSgraGp41+Fpo5nd2rCrRf3plXTiPpvikiURRpQqcD2w2zfSYTddEXixaHW6xPfXEZR2NTxQf2O1/p5InEk0oD6CPhf59yVZlYM4JyrBdwFzI1WcSJVbd23e3nw1QXY+l0AnN67FTdq6rhIXKrMfVAfAxuccwtC23oDRcDgaBQmUpVKSst45Z3lTH53BWWBII0b1ObnF/XiVLVeF4lbkd4HtdI51w0YB3THa/0+Ce8+qMIKv1nEZ7lrdvDgqwvYuHU/AEP6t+eqod2po1XHReJaZe6D2gU8FMVaRKpUQVEJz87K5Y1P1wLQumldbrtEq46LJIqKWr7Pw1vS6CeZWb8qq0ikCszLzeeRyQvZvqeItNQURp/ThbHndaVmDa06LpIofqrlu0hC2b3vIE9MW8xHCzYB0KVtQ267pA8dWjXwuTIRqayKWr7fHctCRI5FMBjkg6838sS0JewrKKZmjTSuHJLF8DM7kaZVx0USUmVavl8C/BfQBTgBuBnIN7N7o1SbSETK33Dbp0tTbrm4Ny0a1/G5MhE5FpGuxTceuBf4B/A/oc15wH3OuXQz+0t0yhM5skAgyBufruHZWbnf3XB73YgenHuybrgVSQaRjqB+CdxkZpOcc78FMLMJzrldeMGlgJKY2rBlHw++uoBla3cCcHqvVtyYoxtuRZJJpAHVCZh/mO0LADXJkZgpLQsw5f2VvDTHKC0LcFy9Wtw0uhenZbfyuzQRqWKRBpQB5wFPlNt+Cd6lPpGoW7lxN/965RvWbN4LwKB+x3PN8B7Uzazpc2UiEg2RBtRvgcnOuZNC3/Nz51xnYBhejyiRqDlYUsZLb+Ux9UOvw23zRpncdnEfendt6ndpIhJFkS519KZzrh/eLL4lwCBgGXCqmX0dxfqkmlu6egcPvvoNm7YdICUFRg7oxBWDs6itDrciSa+ilSSGArPNLABgZkuB8TGqS6q58ssUtW1ej9vH9iGrXSN/CxORmKnoz9DpwA7n3EvAs2a2oIJ9RarM13lbeWjyArbtKiQtNYWLz+3KJed1oUa6likSqU4qCqjWwKWhjzucc0uBZ4EXzWxzLIqT6mVfQTETpi/hvfkbAOjcpgG3j+2rZYpEqqmKljraCvwL+Jdzrj1wGXAlcI9z7n3gOeA1tduQqvDpos38e8oidu87SM30VC6/IItRZ3UiLS3V79JExCeRTpJYC9yDF07dgcvxZvY94px7zcyujl6Jksx27SvisSmL+WSRNyjv3qERt4/tS+umdX2uTET8VumpUGaW65y7F+/eqP/AG1UpoKRSvl/cdTH7CkqoXTON8Rd2Z0j/DqRqcVcRoXKLxTYARuHdnHsusBp4AbgoOqVJstq+u5CHJy9k/rItAPTt2pRbL+5Ds0aZPlcmIvGkwoAKC6WL8VaS2Au8DPyvmc2LfnmSTILBIHO+WMdTM5ZSUFSqxV1FpEIV3Qc1C2+kFABm4K0YMdvMSmNUmySR/B0HeGjSAhau2A7AKT1acNPoXjRukOFzZSISryoaQdUFbgEmmdneGNUjSSYQCDLrkzU8+0YuB4vLqF+nJjeMymZA39YaNYlIhSqaZn5WLAuR5LNp234eePmb71piDOjTmhtysmlQt5bPlYlIItCCZlLlysoCTPtwFS++lUdx6aGWGL05Lbul36WJSAJRQEmVWvvtXh545RtWbtgNwLknt+W6ET3VEkNEKs2XgAqtjD7TzJqFHtcEHsKbiFEG3Gdm9/hRmxydktIAk99dzqvvLqe0LEiThhncdnEfTshq5ndpIpKgYhpQzrkU4Fq8NvHh7gYcXufeBsBs59wmM3sulvXJ0Vm5YTcPvPINa7/15tIM6d+e8Rd2J7N2DZ8rE5FEFusR1N3AhcAfgd+Fbb8KGG9mu4BdoZUqbsRb70/iVHFJGS++lcfUD1YSCELLxnW4bWwfsjs18bs0EUkCsQ6oR83s9865gYc2OOcaAi2B3LD98oDsGNcmlbBszU4eeOUbNm3bT0oKjDqrE+MGZ1G7pt7WFJGqEdPfJkdo03FoVdCCsG0FgNa9iUNFB0t57s1lzPx4NcEgtG1el9vH9lUjQRGpcvHw5+6B0L/hSwpkAvt9qEUqsHDFNh58dQFbdhaQmprCmHM6c+kgR80aaiQoIlXP94Ays13OuXy8SRKbQpuz+OElP/HRgcISnp65lLc+XwdAh1b1uWNsXzq1aehzZSKSzHwPqJCJwF3OuUV4l/x+BTzgb0kCMC83n4cnL2THniLS01K5dFBXRp/ThXQ1EhSRKIuXgPo98A9gKZAKPA486mtF1dzeA8VMmL6Y97/aCEDX4xty+9i+tGtR3+fKRKS68CWgzOwDoGHY4yK8hWlv8aMe+aFPFm3m0dcWsXu/13593NNIstYAAAx+SURBVOBujDyrE2lqJCgiMRQvIyiJA+Xbr/fo2JjbL+lDK7VfFxEfKKDk8O3Xh/VgyGnt1X5dRHyjgKrmduzx2q/Py1X7dRGJLwqoasprv76ep2Ys8dqv107nupE91X5dROKGAqoa2rKzgAdf/ea79uv9urfg5jFqvy4i8UUBVY0car/+3Bu5FBWXUS+zJjfmqP26iMQnBVQ1sWnbfv71yjfkrvHar5/RuxU35vSiYT21XxeR+KSASnKH2q+/8FYeJWq/LiIJRAGVxMq3Xz/npLZcP1Lt10UkMSigklBJaYDJ763g1Xfsu/brt17cmxOzmvtdmohIxBRQSaZ8+/XBp7Xn6mFqvy4iiUcBlSSKS8p4aY4x5YOVBAJBWjTO5PZL+pLdWe3XRSQxKaCSQPn26yMHdOKKwVnUrqXTKyKJS7/BEljRwVImvrmMGeHt1y/pS1Z7tV8XkcSngEpQar8uIslOAZVgCopKeHpmLrM/Wwt47ddvH9uXzmq/LiJJRgGVQOYv28LDkxawfU8R6WkpXDrIqf26iCQtBVQC2FdQzBPTvm+/3qVtQ+4Y25d2LdV+XUSSlwIqzn26aDP/nrKI3fvC2q8P6EiaRk0ikuQUUHFq174iHpu6mE8Weu3Xu3doxO1j+9Ja7ddFpJpQQMWZYDDIh99s4vGpi9lXUOy1X7+wO0P6d1D7dRGpVhRQcaR8+/U+ofbrzdV+XUSqIQVUHAgGg7z95Xqeen0JB0Lt168d0ZPz+qn9uohUXwoon23ZWcBDry5gwYptgNqvi4gcooDySSAQ5I1P1/DsLLVfFxE5HAWUD9R+XUTkpymgYqisLMD0j1bxwuw8iksDNKxXi5su6kX/Xq38Lk1EJO4ooGLkcO3XrxvZk3pqvy4iclgKqCg7XPv1W8b05qRuar8uIlIRBVQUlW+/PuS09oxX+3URkYgooKKguKSMl982Xnv/+/brt13Sh16dm/pdmohIwlBAVbFla3byr1e/YeNWtV8XETkW+q1ZRYoOljJx9jJmzPXar7dpVpc7xqr9uojI0VJAVYFFK7326/k71H5dRKSqKKCOQfn26+1b1ueOsX3p3Fbt10VEjpUC6ijNX7aFhycvZPvuQtLTUhg7yDH67C7USFcjQRGRqqCAqqR9BcVMmL6E9+ZvANR+XUQkWhRQlfDZ4s088lp4+/UsRg7opPbrIiJRoICKwO59B3l06qLv2q/36NiY2y7po/brIiJRpICqwOHar191YXeGqv26iEjUKaCOYMeeQh6ZvIgvc/MB6NOlKbdeovbrIiKxEjcB5ZzrDTwK9AJWA9eY2bxY1xEMBnnny/U8GWq/nlk7nWuG9+T8U9R+XUQkluIioJxzNYHpwP3AAGA0MMc5187M9saqji07C3ho0gIWLFf7dRERv8VFQAEDgRpmdn/o8cvOuVuBscAT0X7xQCDIm5+u4Zmw9us35GRzltqvi4j4Jl4CqjuwrNy2PCA7Fi/+xLTFzPxkDaD26yIi8SJeAqouUFBuWwEQkxkJGbXTadWkDldd2F3t10VE4kS8BNQBoPwbPZnA/li8+M+GdudnQ7vH4qVERCRC8bIEQi7gym3LCm0XEZFqKF5GUO8DKc65/wQewpvF1wuY6mtVIiLim7gYQZlZMTAEL5h2AncCo8xsm6+FiYiIb+JlBIWZLQHO8LsOERGJD3ExghIRESlPASUiInFJASUiInEpbt6DOgZpAPn5+X7XISIilRT2uzut/NeSIaBaAowbN87vOkRE5Oi1BFaFb0iGgJoHnAl8C5T5XIuIiFROGl44/ai9UkowGIx9OSIiIj9BkyRERCQuKaBERCQuKaBERCQuKaBERCQuKaBERCQuKaBERCQuKaBERCQuKaBERCQuJcNKEhFxzvUGHsXr1LsauMbMfnTnsnPueOBJ4FRgK3Cbmb0Ry1qPRiWO7xzgbaAwbPNfzewPMSn0GDnn+gEzzazZEb6ekOfvkAiOL+HOn3NuEPAXoAveOfm7mT12mP0S8txV4vgS7twBOOeGAX8GOuAd399idf6qxQjKOVcTmA68AjQE/gTMcc7VP8zuLwOLgMbA9cDLzrmOsar1aFTy+E4AJplZ3bCPuP4fBMA5l+Kcuw6YA9SsYNeEO39QqeNLqPPnnGsLvAb8Ee9n8zLgHufcBYfZPeHOXSWPL6HOHYBzriUwGfh/ZlYPuBi43zl3wmF2r/LzVy0CChgI1DCz+82sxMxeBpYCY8N3cs51BU4Cfm9mxWb2HvA6cG2sC66kgURwfCEnAgtiWVwVuRu4Ce8XwWEl8PmDCI4vJNHOX3vgRTObamaB0Kj+A+D08J0S+Ny1J4LjC0m0c4eZfQs0NbM3nXOpeOFTCuwL3y9a56+6XOLrDiwrty0PyD7MfuvN7EC5/fpFsbaqEOnxgfdXXFPn3E1ACt6o63dmdjC6JR6zR83s9865gRXsk6jnDyI7Pkiw82dmc4G5hx475xrhLe48sdyuCXnuKnF8kGDn7hAz2+ecywT24GXGX81sRbndonL+qssIqi5QUG5bAZB5lPvFm4jqds6lAxuBqUA34BzgPCCuLzMAmNnmCHZL1PMX0fEl8vkDcM41wPur+gu8S9LhEvbcHVLR8SX6uQOKgDrAycA1zrnyI6OonL/qMoI6AGSU25YJ7D/K/eJNRHWbWSlwbtimlc65PwF/Bf47qhXGRqKev4gk8vkLXQKaDuQC48wsUG6XhD53P3V8iXzuAELHUwzMd849DozEmxBxSFTOX3UZQeUCrty2rND28vsd75zL+In94k1Ex+eca+2cuzc0qeKQmnh/HSWDRD1/EUnU8+ecG4A3qpgGjDGzw9WbsOcukuNL4HN3lnPuq3KbawG7y22LyvmrLiOo94EU59x/Ag8Bo/GmY08N38nMzDm3EPiTc+43QH+8vxROi3G9lRXR8QE7gHFAgXPu//Cmjf4OeCqGtUZNAp+/SCXc+XPOdQJmAnea2YNH2i9Rz12kx0cCnruQBUBr59wvgAeAU/AmPuSE7xSt81ctRlBmVgwMwfvFvRO4ExhlZtucc+Occ+HD0NF414i3AhOAa81sSaxrroxIjy/0l90QYADe/zAfAZOA+3wpvAokw/mrSBKcv1uAenhTr/eHffw1Sc5dRMeXoOcOM9sDDAUuwvvd8jhwnZl9GIvzp466IiISl6rFCEpERBKPAkpEROKSAkpEROKSAkpEROKSAkpEROKSAkpEROJSdblRV6TKOeeeAa6qYJe78Va2fh+oZ2YxWbbHOZcGfAL8zMyWV7BfKvA5cKWZWSxqE6kMjaBEjt4dQMvQx8DQtn5h2+4FPg19fuAw3x8ttwMLKwon+G59tf/Da3QpEnd0o65IFXDO9QQWAx3MbK2PddQG1gPnRHoXv3NuFd5d/x9EszaRytIlPpEoCvV3+u4Sn3MuiNd19Td4C/zOB64A/gu4EtgL/MbMJoa+vx7wD2AMEATeA+6ooD3HpcDu8HByzv0PcAPQFK9v2G/N7M2w75mKNxr8oAoOWaTK6BKfSOz9BfgP4FTgeOBrvGA6GZgCPOacqxva93G8ILsAOAsvpN4K9Rc6nAuB2YceOOdyQq91Bd7q0rOASc65+mHfMxs4r4LnFPGFAkok9h42s/fNbAHeStj78UY1hrd4aAbQwTnXEW9EdLmZzQuNiq7EazM++AjPfRKwNOxxe+AgsC506fH/8Bb+LAnbJxev4VxWlRydSBXRX0wisbcy7PMCYK2ZHXoz+FB/oFpAu9Dn5twP2n1l4o2qZh7muZsD28MeP48303B1qK/P68DTZlYYts+O0L/NKnkcIlGlEZRI7JWUe1y+u+wh6aF9+wJ9wj66Ak8f4XsCQMqhB2a2DTgRb8T1KTAeWBSa1HHIod8DZREfgUgMKKBE4tcyoAZQx8xWmtlK4Fvg73ghdTj5eJMhAHDOXQTcaGZzzOwOvJHXPrweP4c0DftekbihS3wicSrUpfR14Dnn3C3ANuBPeJMr8o7wbV8BvcMepwF/d85twZsxeCrQIvT5Ib2BXfzw0qOI7zSCEolvV+GFyTRgHtAAGGRmu4+w/yy82X4AmNkk4C68Uddy4I/ArWb2Xtj3DABmm5ku8Ulc0Y26IknEOZcJrAUGm9nXEeyfCqzDmyk4N8rliVSKRlAiScTMCvBGS7dE+C0jgdUKJ4lHCiiR5PNPoJcrNze9vNDo6U7g5zGpSqSSdIlPRETikkZQIiISlxRQIiISlxRQIiISlxRQIiISlxRQIiISl/4/QnrtbXEnEGQAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "def plot_velocity(results):\n",
    "    plot(results.v, label='v')\n",
    "    decorate(xlabel='Time (s)',\n",
    "             ylabel='Velocity (m/s)')\n",
    "    \n",
    "plot_velocity(results)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Stopping at 100 kph\n",
    "\n",
    "We'll use an event function to stop the simulation when we reach 100 kph."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "def event_func(state, t, system):\n",
    "    \"\"\"Stops when we get to 100 km/hour.\n",
    "    \n",
    "    state: State object\n",
    "    t: time\n",
    "    system: System object \n",
    "    \n",
    "    returns: difference from 100 km/hour\n",
    "    \"\"\"\n",
    "    x, v = state\n",
    "    \n",
    "    # convert to km/hour\n",
    "    factor = (1 * m/s).to(km/hour)\n",
    "    v = magnitude(v * factor)\n",
    "    \n",
    "    return v - 100     "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>values</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>success</th>\n",
       "      <td>True</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>message</th>\n",
       "      <td>A termination event occurred.</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "success                             True\n",
       "message    A termination event occurred.\n",
       "dtype: object"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "results, details = run_ode_solver(system, slope_func, events=event_func)\n",
    "details"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's what the results look like."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Saving figure to file figs/chap11-fig02.pdf\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjgAAAI4CAYAAABndZP2AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nOzdeXxU9b3/8Vcm+0YgQBJIgBCQL8i+gyyiLe5Ql9b2itati167XWt7be2mtrWttdVr7c+2tnXpYq0tCmLdRTZllR2+LIGwJYSEQELWycz8/jgDhgg4kWTOzOT9fDx4JHPmZOZzHE54+13jAoEAIiIiIrHE43YBIiIiIu1NAUdERERijgKOiIiIxBwFHBEREYk5CW4XcLaMMcnAeKAU8LlcjoiIiIRPPNALWGmtbWz5RNQHHJxws9jtIkRERMQ104AlLQ/EQsApBfjrX/9KXl6e27WIiIhImJSVlTFnzhwIZoGWYiHg+ADy8vIoKChwuxYREREJvw8NUdEgYxEREYk5CjgiIiIScxRwREREJOYo4IiIiEjMiYVBxh+purqa8vJyvF6v26W0i/T0dAoKCvB4lE9FREROJeYDTnV1NQcPHiQ/P5/U1FTi4uLcLums+P1+9u/fT0VFBTk5OW6XIyIiEpFivgmgvLyc/Px80tLSoj7cAHg8HnJzczl69KjbpYiIiHykQCDgyvvGfMDxer2kpqa6XUa7SkxMpLm52e0yRERETsvb7OPhZ9dw/Q9fobSiNuzvH/NdVEBMtNy0FGvXIyIiseVYvZcHnlzB+h0VpCbHkxAf/vaUThFwREREJDwqjtTzoz+8S0lZDd0yk/nhFybRs1v4e1IUcERERKRdlJRW86M/vEvF0QYKcjL40Rcnk5ud5kotCjgiIiJy1tbvOMRP/7yC2oZmhhRm8/1bJ5KZluRaPTE/yDha/Pvf/2bs2LEcPHgQcHZHnzFjhmZLiYhIxHtnzT5++Pt3qW1oZvLwXtx/23muhhvohC049z7xHqu2HAzb+40bkssPvzDpI8+7+uqrefPNN7n33nv51re+xUMPPcTjjz9OVlZWGKoUERFpu0AgwNyFO/jzS5sBmDWtiFtnDyPe4/5kmE4XcCLZfffdxxVXXMHNN9/MDTfcwIQJE9wuSURE5JR8Pj+/f2EDLy/bDcDNVwzlqhkDImamb6cLOKG0prile/fuXHTRRTz77LPMnj3b7XJEREROqaGxmV/8ZRUrNx8kId7Dnf81hmmj890u6yQagxNBNmzYwLx587j88sv53ve+h9/vd7skERGRk1RVN/Cd/7eUlZsPkpGayI9vOy/iwg0o4ESMhoYGvv3tb3Pbbbfx05/+lMrKSp544gm3yxIRETlh78Ea7np0MTv2HiE3O40HvzaNoUXd3S7rlBRwIsRDDz1ESkoKt956KykpKdx///08+uijbN261e3SRERE2Lizgm89upjyw3UM6tuVB782jYKcTLfLOq1ONwYnUt1zzz0nPZ44cSIbNmxwqRoREZEPvLNmHw8/+z7NPj8Th+Zx1/VjSUmK7AgR2dWJiIiIawKBAM+9uY2//MfpTbhiSn++cOXwiJgG/lEUcERERORDvM1+Hnt+LW+u3EtcHNw6exizpxVFzDTwj6KAIyIiIic5VtfEA0+tZP2OCpKT4rlrzlgmDevldlltooAjIiIiJ5RV1nLvE++xr/wY3TKT+f6tEzmnTze3y2qzThFwAoFA1DSphSIQCLhdgoiIxKCtJYf58Z+Wc/RYE33zMvnhrZPIcWk38LMV8wEnMTGR+vp60tKi8wM6Fa/XS0JCzH90IiISRkvXHeBXf1tNU7OfUYN6cvfnx5Oemuh2WR9bzP8rmZOTw/79+8nPzyc1NTXqW3L8fj8HDx7UJpwiItIuAoEAz7+1nadf3gLARRP7cfs1I0iIj+6l8mI+4HTp0gWAAwcO4PV6Xa6mfaSnp9OjRw+3yxARkSjnbfbzm3+u5a1Vzkypmy4/l6tmDIz6xgDoBAEHnJBzPOiIiIgIVNc28dMnV7CpuJKkxHjumjOGycN7u11WuwlrwDHGzAR+BpwDlAMPWmt/Z4xJBmqAphanL7PWXhTO+kRERDqDfeU13PfH5ZRW1JLdJZnv3zKJgX26ul1WuwpbwDHG9AH+BdwIvAiMBV41xuwGKoHD1tq8cNUjIiLSGa3fcYgHnlzJsXovRb2z+P6tE+nRNdXtstpdOFtwCoG/WWvnBh+vNMYsBKYA+4G1YaxFRESk03lteQm/fX4dPn+AiUPz+OacsaQmx+ZolbBdlbV2MbD4+GNjTDYwDXgGuATIMcasB3KBRcA3rLX7w1WfiIhIrPL5Azz50iZeeGcnAFeeP4CbrhgaFXtKfVyuzAEzxmQB84DlON1VtcBS4BOAAeqBuad9AREREQlJXYOXH/9pOS+8s5N4TxxfvXYUt84eFtPhBlyYRWWMGYQTajYDc6y1fuDOVufcCRwyxvSx1u4Nd40iIiKxoKyylvv/tJw9ZTVkpiXynZsmMHxA51hmJKwtOMaY6TitNi8An7bWNgSP32eMGdLi1KTg14Zw1iciIhIrNhVX8s1HFrGnrIaCnAwe+vr5nSbcQHhnUQ0AXgLusdY+2urpEcA4Y8x1wcePAAustYfCVZ+IiEiseGPFHh57fi3NvgBjTA7fvmFcVG+78HGEs4vqDiATeMAY80CL448BtwL/B+wI1rQA+FIYaxMREYl6Pn+ApxZsZu7CHQDMnlbELbOGEh/l2y58HOGcRXUnrcbatDInXLWIiIjEmtp6Lw/+ZRWrt5YT74njy1eP4NLJhW6X5ZrYnPwuIiLSiRw4dIz7/7ScfeXHyExL4js3jmf4wM4z3uZUFHBERESi2Npt5fz86VUcq/fSLy+T790ykbzu6W6X5ToFHBERkSgUCASYv6SYP87bhD+4MvGd140hLaVzDSY+HQUcERGRKONt9vO7uet59b0SAD7ziXO4/pIheGJ88b62UMARERGJIlU1DfzsqZVs3nWYpAQPX/3saGaMKXC7rIijgCMiIhIltu+t4qd/XkHF0Qa6Z6Xw3ZsmMKhvN7fLikgKOCIiIlFg4eq9PPrcWpqa/QwpzOY7N46nW5cUt8uKWAo4IiIiEaz14n0XTezHbVcPJzEh3uXKIpsCjoiISIQ6VtfEg39ZzRrrLN73xSuHc9l5hcTFaTDxR1HAERERiUB7yqr58Z9XUFpRS5f0JO6+cXyn2izzbCngiIiIRJhl6w/w8LNrqG/0UdQ7i3tunkBOdprbZUUVBRwREZEI4fMH+OsrW/jnm9sBmD46n69+ZhQpyfrnuq30X0xERCQCHKtr4sG/rmbN1nI8cXDzrKF8avoAjbf5mBRwREREXFZSWs1P/ryC0spaMtOS+N/Pj2PkOT3dLiuqKeCIiIi4aMm6/Tzy7Ps0NPkoys/inps03qY9KOCIiIi4wOfz8/TLW/h3cH2bGWMKuOMzI0lJ0j/N7UH/FUVERMLsSE0jD/5lFet3VODxxHHLrKHMnlak8TbtSAFHREQkjLbtqeKBJ539pLpmJvO/N4xjmNa3aXcKOCIiImEQCAR49b0Sfjd3A80+P4P7dePuG8fTPSvV7dJikgKOiIhIB2vy+nj83+t5fcUeAC6f0p9bZw8jMcHjcmWxSwFHRESkAx08XMfPnlrBjn1HSUrwcMdnRnLhuL5ulxXzFHBEREQ6yKotB/nV31ZTU+clNzuN7940gaL8LLfL6hQUcERERNqZzx/g2dcs/3jDEgjAuCG53HndGDLTktwurdNQwBEREWlHR4818tBfV/P+tkPExcH1lw7mMxcOwuPRFPBwUsARERFpJ9v2VPHAUyupOFJPl/QkvnX9WEYNynG7rE5JAUdEROQsBQIBXl62myde3ECzL4Dp1427Pz+eHl01BdwtCjgiIiJnoa7By2PPr2PR+/sBuGJqf26ZpSngblPAERER+ZhKSqt54KmV7D90jJSkeL567Simjy5wuywhzAHHGDMT+BlwDlAOPGit/Z0xJgn4DfBpwAf8ylr7QDhrExERaYs3V+7ht/9aT5PXR7+8TP738+Ppk5vpdlkSFLaAY4zpA/wLuBF4ERgLvGqM2Q3MAAwwAMgCXjHG7LfWPh2u+kRERELR6PXxuxarEl84rg+3XzNCu4BHmHB+GoXA36y1c4OPVxpjFgJTcELPTdbaKqDKGPNL4MuAAo6IiESM/YeO8bOnVrK7tJqkBA+3XT2CmRP7uV2WnELYAo61djGw+PhjY0w2MA14BugFbG5x+lZgeLhqExER+SiL3t/Hb/65jvrGZnr3SOfuG8fTv7dWJY5UrrSnGWOygHnAcmB18HBdi1PqgLRw1yUiItJao9fHEy9u5JV3dwMwZWRvvnbtKNJSEl2tS84s7AHHGDMIZwzOZmAOcHyRgJaLBaQBx8JcmoiIyEn2ldfw86dXsbu0msQED1/41DAunVxIXJxWJY504Z5FNR0n3DwOfNdaGwAajDFlOIOM9wdPHczJXVYiIiJhtXD1Xh57fh0NTT5690jnfz8/XhtlRpFwzqIaALwE3GOtfbTV088APzTGrAcygLuAR8JVm4iIyHENTc384YWNvLa8BIDpo/O549Mj1SUVZcLZgnMHkAk8YIxpucbNY8APgIeATYAH+D1OK4+IiEjY7Cmr5hfPrKKkrIakBA9fumo4F03spy6pKBTOWVR3Anee4ZQ7gn9ERETCKhAI8NryPfz+hQ00eX3k93S6pDRLKnqFHHCMMbk4i/Pl4Kw2XAassdZWdlBtIiIiHa623tlLavFaZxjoheP6cNvVI0hN1sJ90eyMn54xJgG4DvgGMBJoAqqAeCA7eM5y4LfAs9Zaf4dWKyIi0o627aniwb+soqyyjpSkeP770yO5YGwft8uSdnDarU6NMecD64HPA38EBgFp1tre1tpcIAkYDfwN+Aqw1Rgzo8MrFhEROUt+f4C5C3fw7UcXU1ZZR1F+Fg/fOUPhJoacqQXnm8BnrbUbTvVkcIr3xuCf3xpjRgP3AQvbu0gREZH2cqSmkYefXcPqreUAzJ5WxE1XnEtiQrzLlUl7Om3AsdbObssLWWvfB2addUUiIiIdZM3Wcn797BqO1DSSmZbI1z87monDerldlnSAtgwyTgP6A8mtn7PWrmnPokRERNqTt9nH0y9v4YV3dgIwfEAP7rxuDD26pn7ET0q0CingGGOux1mXJhVovRhAAGfQsYiISMTZV17Dg39ZTfH+o3g8cVx/yWCuvuAc4j1a2yaWhdqC8wDOQONfAQ0dV46IiEj7CAQCvL7CWdumsclHXvc07pozFtMv2+3SJAxCDThdgN9Ya0s6shgREZH2UFPXxGPPr2PpugMAzBhbwO1Xj9B2C51IqAHnGeAm4J6OK0VEROTsrd9xiF/9bQ2VRxtITY7n9mu0tk1nFGrAeRBYY4yZA+wGTlrQz1p7YTvXJSIi0ibeZh9/+c9W5r6zg0AABvfrxp3XjaVXj3S3SxMXtKUF5xiwAKjruHJERETabu/BGn75l9UUH3AGEv/XzEFc+8lBxMefdj1biXGhBpzxwERr7fqOLEZERKQtAoEALy/bzZ/mbaSp2U9e9zS+ed1YBhdqIHFnF2rAsUDXjixERESkLaqqG/i/59ayastBAD4xvg9funK4BhIL0LZp4k8aY34D7AS8LZ+01r7c3oWJiIiczrsbSvnNP9dSXdtEemoiX/nMSKaOzHe7LIkgoQacvwe//vIUz2mhPxERCYu6Bi9/eGEjb6zcA8Coc3ry9c+N1orE8iEhBRxrrUZpiYiIqzYVV/Krv6+h/HAdSQkebrpiKJdP6Y9HKxLLKZw2uBhjzm/rixljNF1cRETalbfZx5MvbeI7v11C+eE6BhRk8fCdM5g1rUjhRk7rTC04/2OMuRv4P+ANa633VCcZYxKAK4Cv4Ewhf6vdqxQRkU5p14Gj/Prva9h1oBpPHFz7yUF8dqYhMUEdC3Jmpw041torjTFXAT8D+hljFgKbgAqcDTd7AiOBycAe4H5r7fMdXrGIiMQ8n8/Pvxfu4G+vbqXZF6BX93TuvG6Mpn9LyM44BsdaOxeYa4yZAVyGE2ZycVYyLgNWAw9Yaxd3cJ0iItJJ7D90jF//fQ22pAqAy84r5KYrhpKaHOq8GJHQBxkvBBZ2aCUiItKp+f0BFizdxZMLNtPk9dE9K4WvfXY0Y0yO26VJFFIcFhER15VX1fHIs++zfkcFABeMLeBLVw4nIy3J5cokWingiIiIawKBAK8t38Mf522kvrGZrIwk/vuakZw3orfbpUmUU8ARERFXVByp59Hn1rLGlgMwaVged3x6FF0zk12uTGKBAo6IiIRVIBDgjRV7eGLeRuoamslMS+RLV43g/NH5xMVpXRtpHyEHHGNMDjACSMSZJn6C9qISEZFQVB51Wm1Wb3VabSYOzeOOT4+kW5cUlyuTWBNSwDHG3Ar8FifctNbmvaiMMROAl6y1OcHHyUAN0NTitGXW2ova8roiIhKZAoEAb67cyxMvbqC2oZmM1ES+dNVwZowpUKuNdIhQW3C+BfwB+I61tubjvpkxJg64lQ9v2jkcOGytzfu4ry0iIpGpvKqOx/657sRYmwnn5nHHZ0aSrVYb6UChBpw+wCNnE26C7gUuB34MfK/F8bHA2rN8bRERiSB+f4BXl5fw5/mbqG90Wm2+eOVwLhirVhvpeKEGnNeATwDbz/L9HrfW/iC4MnJLY4AcY8x6nJWSFwHfsNbuP8v3ExERF5RV1vLoc2tPrGszeXgvbr96hMbaSNiEGnDWAb8yxswGtnHyWBmstd8O5UWstQdO81QtsBS4D/DibPA5F5gQYn0iIhIBfP4AC5YU8/R/ttDY5CMrI4nbrh7BlBG91WojYRVqwDkfWA6k4myw2VLgbIuw1t7Z8rEx5k7gkDGmj7V279m+voiIdLw9ZdU8+txatgb3kDp/dAFfvHIYWRla10bCL9S9qC7oyCKMMfcBf7fWbgkeOr42d0NHvq+IiJw9b7Of59/aznNvbKPZ5ye7SzL/fc1IJg7r5XZp0om1ZR2cXOArwFDAA2wB/mCtLW6HOkYA44wx1wUfPwIssNYeaofXFhGRDrK15DCPPreWPWXOHJSLJ/XjpiuGkpF6qlVFRMIn1HVwJgCvA3uBZTgL/V0BfNUYM8Nau+os67gVZ9zNjmBNC4AvneVriohIB6lvbOYv/9nC/CXFBALQq0c6X/3MKIYP7OF2aSJA6C04DwF/B2631p4Yc2OM+Q3wINCmLixr7UKga4vHlcCctryGiIi4Y/XWg/z2+XWUV9Xj8cRx9YwB/NfFg0lObNOaryIdKtSAMw74QstwE/QosLp9SxIRkUhUVdPAEy9sZNFaZwWPAQVZfPUzoxhQ0PUjflIk/EINOKVAIWBbHS/C2WJBRERilN8f4PUVe/jzS5uorfeSlBjPnIsNn5o+gPh4j9vliZxSqAHnGeD3xphvAO8Fj00Gfh18TkREYtDegzU89vw6NhVXAjBmcA63Xz2CvO7pLlcmcmahBpyfAL2B53BmUMXhLMj3KHBPx5QmIiJu8Tb7eP7N7Tz35naafX66ZiTzhU8NY/rofC3YJ1Eh1HVwmoAvGmPuAgxQD+yw1tZ3ZHEiIhJ+67Yf4v/9az37Dx0DYOaEvtw8ayiZaUkf8ZMikeO0AccYcxnwurXWG/y+tT7GGACstS93UH0iIhImR2oa+eP8jSxcvQ+A/J4Z3PGZkQwfoKnfEn3O1ILzEpAHlAe/P50AoLmBIiJRyu8P8NryEp5csNkZRJzg4dqZg7h6xkASE/TrXaLTaQOOtdZzqu9FRCR27DpwlMeeX4cN7h81xuRw29Uj6NVDg4gluoW6kvFbwNXW2iOtjvcEXrHWju2I4kREpGPUNXj526uW+UuK8fsDZHdJ5otXDteu3xIzzjQGZwZwbvDh+cCXjTGt17wZAgzomNJERKS9BQIBFq/dzx/nbeRwdSOeOLhian9uuHQIaSnaP0pix5lacCqBu3CmhMcBdwC+Fs8HgGPANzusOhERaTd7D9bw+L/Xs35HBQCmbzduu2YEA7USscSgM43B2YCzUjHGmLdxuqiqwlWYiIi0j4bGZp593fLiop00+wJkpiVx4+XnMnNCXzwedUdJbDpTF1WatbYu+PDy48dOdW6L80REJEIEAgGWrS/liXkbqThST1wcXDypH5+/7Fy6pGtNG4ltZ+qiqjHG9LLWluN0RbXeaBOcritNExcRiTB7yqr5/QsbWLfd6Y4aWJDF7deMZFDfbi5XJhIeZwo4FwKHg99fEIZaRETkLNXWe/n7a5aXlhTj8wfITEvkhkuHcNGkQuLVHSWdyJnG4Lxzqu8BjDFJwAhgm7W2uuPKExGRUPj9Ad5evZcnF2zmSI0zO+rS8wq5/pIh6o6STinUdXAGAn8E/hdYDyzDCThHjTGXWmvfO9PPi4hIx9m+t4rfzd1wYrG+IYXZfPmq4QzQ7CjpxELdTfxRoAbYDdwAFOBsunkz8CvgvI4oTkRETq+quoGnX97CGyv3ANAtM5mbZw1lxpgCLdYnnV6oAWcaMNpaW2aMuRJYYK3dboz5A/CNjitPRERa8zb7mLeomH+8Yalv9JEQH8enpg/g2k8O0mJ9IkGhBpwGINEYk46zqvEtweN5wNGOKExERE4WCARYvqmMP83bRGllLQATh+Zxy6yh9O6Z4XJ1IpEl1IDzKvAHnG6qOmC+MeYTwCPAvA6qTUREgkpKq3li3kbWbjsEQJ/cDL7wqeGMMTkuVyYSmUINOF8GfgL0Ay631tYaY8YDC4FvdVBtIiKdXlVNA3971fLae7vxByA9NZE5Fw/m0vMKSYj3uF2eSMQKKeBYa48BXwcwxnQxxnS11v6sQysTEenEmrw+Xly0k3++uZ36xmY8njiuOK+Qz11kyMpIdrs8kYgXagsOxpjbge8CvYOPy4FHFHRERNpPIBBgyboDPLlgM+WHnV1wxg3J5ZZZQ+mTm+lydSLRI9R1cO4Cvo/TTbUEZ4uGKcDdxph6a+0jHVeiiEjnsLXkMH+at4ktu51F5PvlZXLL7GEaZyPyMYTagnMHcJu19u8tji01xpQAP8YZbCwiIh9DaUUtT728maXrDgCQlZHEnEuGcNGEvsRrnI3IxxJqwOkJrDzF8dU4i/6JiEgbVdc28Y/XLS8v20WzL0BSgodPnT+Aay44h/RUrWcjcjZCDTgbgc8AD7Q6/llga1vf1BgzAXjJWpsTfJwE/Ab4NOADfmWtbf1eIiIxocnr46UlxTz3xjZqG5qJi4MLx/Xh+kuG0LNbqtvlicSEUAPOD4AFxpjJwLvBY5OBS4CrQ30zY0wccCvwy1ZP3Yuz9cMAIAt4xRiz31r7dKivLSIS6Xz+AO+s2cdfX9lCeVU9AKPO6cnNs4ZSlJ/lcnUisSXUaeKvBRf2+yrOXlT1wBZgvLV2XRve717gcpxxO99rcfxG4CZrbRVQZYz5Jc7aOwo4IhL1AoEAq7eW89SCzewurQacAcQ3zxrKGJOjfaNEOkDI08SttYuARWf5fo9ba39gjJlx/IAxpivQC9jc4rytwPCzfC8REddt21PFky9tZsPOCgB6dE1lzsWDuWBcH+I9CjYiHeW0AccYkwY8jDMuphGYC9xtra3+uG9mrT1wisPHN1Cpa3GsDkj7uO8jIuK2/YeO8czLW1i63vm1l5GayLWfHMTlU/qTlBjvcnUise9MLTj3ArOAX+AM/P0K0B1nYHF7qg1+bTmyLg041s7vIyLS4Q5V1fOPNyyvr9iD3+/MjJo9fQDXXHgOGZoZJRI2Zwo4nwaus9a+DWCMeQdYZIxJtNZ626sAa22VMaYMZ5Dx/uDhwZzcZSUiEtGOHmvkn29u5+Vlu/A2+/F44rhoYj+uu9jQPUszo0TC7UwBp4CTp4CvBDxALrCvnet4BvihMWY9TpfVXWjxQBGJArX1Xl54ZycvLtpBfaMPgGmj8rnuYkNBjrZWEHHLmQJOPE7XFADW2oAxphFI6oA6fgA8BGzCCVG/Bx7vgPcREWkXDU3NvLx0N8+/tY2aOqdRe9yQXG64dIimfItEgJBnUbUna+1CoGuLxw0420Hc4UY9IiKhavL6eOW93Tz/5naqahoBGFrUnRsuHcLQou4uVycix31UwLnJGNNysG8CcL0xpqLlSdba37Z7ZSIiEcTb7OeNFSX8441tVB5tAGBgn65cf8lgrWUjEoHOFHD2ALe3OlYG3NzqWABQwBGRmNTs8/PWqr3843V7YvXhwl5dmHPJYCYOzVOwEYlQpw041trCMNYhIhJRfD4/C9fs4x9vbKO0wlnNok9uBtddPJjzhvfGo0X6RCKaK2NwREQilc/n5+3V+3jujW2UVjrBpnePdP7rIsO00QVafVgkSijgiIjgdEUtXL2Xf7yxjbJKZ2H13j3S+exMw/mj84mP97hcoYi0hQKOiHRqzT4/b6/ay3NvfhBs8ns6wWb6KAUbkWilgCMinVKT18frK/bwr7e3cyg4eDi/Zzqfm6muKJFYoIAjIp1KQ2Mzr7y3m7kLd3C42lnHpk9uBtd+YpCCjUgMUcARkU6hrsHLgqW7eOGdnVTXNgFQ1DuLa2cOYvKwXpoVJRJjFHBEJKYdqWlk3uKdvLx0F7UNzQCYvt24duYgxg/J1To2IjFKAUdEYtLBw3XMXbiD15eX0NTsB5wtFT43cxAjz+mpYCMS4xRwRCSmlJRV8/xb21n0/n78/gAAE87N49MXnsOQ/tkuVyci4aKAIyJRLxAIsHnXYeYu3MHyTWUAeDxxzBhbwDUXnENhry4uVygi4aaAIyJRy+cPsGJTKf96ewe2pAqAxAQPMyf05aoZA8nrnu5yhSLiFgUcEYk6jV4fb63aywsLd3AguE9URmoil0/pz+VT+9MtM8XlCkXEbQo4IhI1jh5r5JV3d/PSkl0cOeasYZOTncaV0wcwc0JfUpL1K01EHPptICIRb+/BGl5ctJO3V+09MSNqQEEWV88YyJQRvbWdgoh8iAKOiESkQCDAuu2HeHFRMau2HDxxfNyQXK48fwAjBvbQVG8ROS0FHBGJKE1eH4ve38eLi6OFXxIAACAASURBVIrZXVoNQFJiPJ8Y14dZ04rok5vpcoUiEg0UcEQkIlQereflZbt55d3dJ7ZS6JaZzOVT+3PJpEKyMpLdLVBEoooCjoi4JhAIsHV3FfOXFLNs/QF8wYX5ivKzmD2tiOmj80lMiHe5ShGJRgo4IhJ2TV4fS9btZ/7iYnbsOwo4C/NNHdmbWdOKGFKYrfE1InJWFHBEJGzKKmt55d3dvLZ8DzV1TjdUZloSl0zux6WT+9OzW6q7BYpIzFDAEZEO5fcHeH9bOQuW7mLVloMEnF4oivKzuHxKf84fU0ByorqhRKR9KeCISIc4eqyRN1fu5ZV3d1Na6aw2nBDvYeqo3lx+Xn9Mv27qhhKRDqOAIyLtJhAIsKm4klfeLWHp+gM0+5xF+Xp2S+XSyYXMnNCPrpmaDSUiHU8BR0TO2rG6Jt5atZdX3tvN3oPHAIiLg7GDc7hkciHjz80j3qPWGhEJn4gJOMaYW4DfAY0tDt9hrX3KpZJE5AwCgQCbdx3mteUlLFm7/8QWCt0yk5k5sR8XTexHbnaay1WKSGcVMQEHGAM8ZK292+1CROT0qmoaeGvlXl5fUcL+Q7Unjo8a1JNLJhcycWgeCdobSkRcFkkBZyzwiNtFiMiH+fwB3rflvLa8hBWbyk4syJfdJZlPjO/LzAn96NUj3eUqRUQ+EBEBxxgTD4wAbjDG/AqoA54Afm6tDbhanEgntvdgDW+t2stbq/ZyuLoBcBbkmzg0j4sm9WOsydFO3iISkSIi4AA9gVXAU8DVwBDgRaAa+K2LdYl0OsfqvSxZu583Vu7BllSdON6rRzozJ/TlE+P7kt0lxcUKRUQ+WkQEHGttGXB+i0NrjTGPAteggCPS4Xw+P+u2V/Dmqj28t6H0xIDh1OQEpo7szScn9NX2CSISVSIi4BhjhgLXWmt/2OJwEtDgUkkiMS8QCFC8/yhvr97Hovf3UVXzwQTGEQN78MkJfZk8rBcpyRHxa0JEpE0i5TfXEeCbxph9wB+B0cDXgK+4WpVIDDpUVc877+/j7dV72VNWc+J47x7pXDCuDxeO7UOOpneLSJSLiIBjrd1vjJkN/AL4NVAB3G+tfd7dykRiQ3VtE8vWH+Cd9/exqbjyxH5QmWlJTB+dzwVjCxjUV1sniEjsiIiAA2CtfQsY53YdIrGivrGZ5ZvKeGfNPt635SemdicmeJhwbh4XjC1gzOBcEhM0C0pEYk/EBBwROXtNXh+rt5azZO1+lm8uo7HJB4AnDkYP6sn00QVMHt6L9NRElysVEelYCjgiUa7J62ONLWfJ2gOs2FxGfWPzieeGFGYzfXQ+U0b2plumpnaLSOehgCMShZq8Pt635SxZf4DlG08ONQMLspgyMp/po/I1WFhEOi0FHJEoUd/YzKotB3l3QymrtpRR3+g78VxRfhZTR/Zm6sh8bZkgIoICjkhEq6lrYsWmMt7dUMoaW443uAAfQFHvLKaO6s2Ukb3p3SPDxSpFRCKPAo5IhCk/XMd7m0pZvrGMTcWVJ2Y/gTOmZvLwXkwe3ou87mqpERE5HQUcEZcFAgF27jvK8k1lLN9Uyq4D1See83jiGHlODyYP782kYXl0z0p1sVIRkeihgCPigoamZtbvqGDV5oOs3FxGxdEPdiVJTY5njMll4rA8xg3JJTMtycVKRUSikwKOSJiUH65j5ZaDrNpykPXbD53Y0BIgu0sKE4fmMXFYHiMG9iAxId7FSkVEop8CjkgH8Tb72Fx8mNW2nDVbD1LSYt8ngHP6dGX8kFzGnZvLgPyueDzaJkFEpL0o4Ii0owMVx3h/azmrbTnrd1ScWEkYIDU5gdGmJ+OH5DJ2cC7dumjhPRGRjqKAI3IWauqaWL+jgnXbDrF22yFKK2tPer6wVxfGDs5htMnh3P7dte+TiEiYKOCItEGT18eW3YdZt90JNDv2HTmxMzdARmoio00OY0xPRpsczXoSEXGJAo7IGXib/WzfW8WGHRWs31HB1t2HTxocnBAfx5DC7owc1INR5/RkYJ9uxGssjYiI6xRwRFrwNvvZue8IG3Y6gWbL7sMnjaMBp9tp1KCejB6Uw7n9s0lJ1m0kIhJp9JtZOrWGxma2lhxmU/FhNu+qZGtJFU3ekwNNn9xMRgzswfCBPRhW1J2sjGSXqhURkVAp4Eincri6ga27D7NltxNodu47etJWCAB9cjMYWtSDEQN6MGxgd7plaraTiEi0UcCRmOXz+dldWh0MNFVsKTlM+eG6k87xxMHAgiyGFvVgaFE25/ZXC42ISCxQwJGYEAgEOFRVz7a9VWzbc4Rte6rYse/Ih8bPpCbHY/pmM7gwmyGF2Qwu7EZaSqJLVYuISEdRwJGoVFXTwM59R9mx7wjb9xxh294qjtQ0fui8vO5pDDkRZrLpm9dFs5xERDoBBRyJaIFAgMPVDezcf5Sde4+wIxhqDlc3fOjczLREBvXtduLPOX26qrtJRKSTUsCRiOFt9rOvvIZdB46y60D1ia/VtU0fOjc1OYEBBVkMyO/KwD5dMX27kdc9jbg4tc6IiIgCjrjA7w9QXlXHnrIaSsqqKSl1vu4rr6HZF/jQ+empiRT1zmJAQRYDC5xA06t7ujanFBGR01LAkQ7T7PNTVlnLvvJjwT817CmrYe/BGhpaDf49rlePdPr37kL/3lkU9c6isHcXenZNVcuMiIi0iQKOnJXjY2RKK2o5UFHLgUPHOFBRy77yGkorak/ZIgOQ3SWZvnld6JuXSb+8LvTLy6RPbqZmNImISLtQwJGP5G32c6iqjrLDdRw8XMfBylrKKus4UHGM0ora07bGAOR0S6UgJ5OC3AwKcjLpk5NB37wudElPCuMViIhIZ6OAI9Q3NnOoqo5DR+o5VFVPxZF6Dh2pd8LM4Toqj9aftGN2a5lpifTukUGvHun07pFOr54Z9MnJIL9nhvZpEhERV0TMvz7GmJHA48AIoBi4xVq70t2qoluT18eRY40cqWmk8mgDh6uDf442UHm0nsPVDVQebeBYvfeMr+OJgx7dUsnNTiMvO53c7mnkZqfRu0c6vXtmkJmm1hgREYksERFwjDFJwIvAw8B04BrgNWNMP2tttavFRQi/P0B9YzM1dU1U1zZRU9dETa3zfXXw2NFgmDlS08iRY43UNTSH9NqJCR56dk2lZ7dUenZNC351Hudmp9OjayqJCZ4OvkIREZH2ExEBB5gBJFprHw4+ftYY8xXgs8AfXKvqY/D5A/h8fpp9fpq8fpqafTR5fR/6vqGpmfqGZuqbmmlo9J14XNfYTF2Dl9r6ZmrrvdQ2eKmt91LX4MV/hm6iU4n3xJGVkUzXzGS6Z6WQ3SWF7l1SyM5KoXtWKtldnGNZGUmapSQiIjElUgLOucCWVse2AsNdqOWEugYv9/9pOYeq6gkEAgSAgD+AP+DMHgoEwOd3wkyzzwk2bQ0hbZGSFE+X9CQy05PITEuiS1rSSY+7ZibTNRhoumYmk5GaqOAiIiKdUqQEnAygrtWxOiDNhVpOqKnzYkuq8Db72/RzCfEeEuLjSEqMJynB43xNjCc5MZ7EROdxalICKcnxpCYnkJqcQMrxx0kJpKcmfvAn5fjXBOLj1U0kIiISikgJOLVAaqtjacAxF2o5ITc7jad/dAnH6pytAjxxccTFxeEJ5oy4uDjiPXEkJniIj/eQ4InD44lTq4mIiIjLIiXgbAb+p9WxwcDTLtRykozURDJStficiIhINImUgPM2EGeM+R/gNzizqEYAc12tSkRERKJSRAzqsNY2AZfiBJvDwD3AldbaQ64WJiIiIlEpUlpwsNZuBKa6XYeIiIhEv4howRERERFpTwo4IiIiEnMUcERERCTmKOCIiIhIzImYQcZnIR6grKzM7TpEREQkjFr82x/f+rlYCDi9AObMmeN2HSIiIuKOXsDOlgdiIeCsBKYBpYDP5VpEREQkfOJxws3K1k/EBQIduP21iIiIiAs0yFhERERijgKOiIiIxBwFHBEREYk5CjgiIiIScxRwREREJOYo4IiIiEjMUcARERGRmBMLC/21mTFmJPA4MAIoBm6x1n5okSBjTF/gj8AkoBz4qrX25eBzccD9wJeAJODPwLestc1huYgQtOE6xwK/Dp5XDTwB3G+tDQSf3wN0B44vmrTfWms6/gpC04brvBB4Hahvcfjn1tr7Y+XzNMZMA/7T6keTgV3W2kHBcyL68zzOGDMBeMlam3Oa56P6/jwuhOuM6vvzuBCuM6rvz+POdJ3Rfn8aY2YCPwPOwbnnHrTW/u4U50XEvdnpAo4xJgl4EXgYmA5cA7xmjOlnra1udfqzwLvA5cBU4AVjzChrbTHOh3M1MAZoBOYC3wXuC8uFfIRQr9MYkwYsAH4MXAAUAa8CZcDvjTE9gHygi7W2NrxX8dHa+HmOAf5prf3cKV4qJj5Pa+1iIKPFz/XBWeHzK8HHEf15wolfgLcCv/yIU6P2/oTQrjPa709o0+cZtfcnhHad0Xx/Bmv9F3Ajzu+iscCrxpjd1tpXW50eEfdmZ+yimgEkWmsfttZ6rbXPApuAz7Y8yRgzCBgH/MBa22StfQuYh/MXGJwP+WFr7T5r7SHgR8CXw3QNoZhBCNcJ9AHetdb+xlrrs9ZuB17A+UsJzl/i7ZF2s7Uwg9CuE5xrWXua14mVz7O1PwHPWGtfCz6O9M8T4F7gdpx/1E8pBu5PCOE6if77E0K7Toju+xNCv86Woun+LAT+Zq2da631B1uPFwJTWp4USfdmZww45wJbWh3bCgw/xXl7Wv1Fa3neucDmVs/1NsZkt2OtZyOk67SOq44/DrYUXAq8Hzw0BvAYY1YYYw4ZY141xgzpwLrbKtTPE5xr+aQxpsQYs8cY86AxJrnF60T959mSMebK4M/9oMXhSP88AR631o4FVp3hnGi/PyGE64yB+xNC+zwhuu9PCP06gei7P621i621tx1/HPxvP40P/i4eFzH3ZmcMOBlAXatjdUBaG89r/fzx71u/jltCvc4Tgr9M/h487/HgYR+wAqdJsR/OX+b/BJvOI0FI12mMSQD24TSHDgEuBD6J0xd8qteJ+s8TuAf4mbW25ZiGSP88sdYeCOG0aL8/Q73OE6L0/gzpOmPg/mzz50mU3p8AxpgsnFaZ5TjdVS1FzL3Z6cbgALVAaqtjacCxNp7X+vnjH07r13FLqNcJgDEmD6d/1Q988vhNZ639RavzvgP8N05T6uJ2rvnjCOk6gwPYPtHi0A5jzE+AnwPfPsXrRPvnOQIYCjzV8ngUfJ6hivb7s02i+P4MSQzcn20SzfdnsAvqRZxWmDnWWn+rUyLm3uyMLTibgdYj0gdzcpPZ8fP6GmNST3Ne69cZDJRaa4+0Y61nI9TrxBhzLs5Atx04vzyrWjz3DWPM1Banx+ME44Z2r/jjCek6jTH5xphfBpv4j0vig+uImc8z6FPAf1oPtI6CzzNU0X5/hizK78+QxMD92VZReX8aY6bjtNq8AHzaWnuquiLm3uyMLThvA3HGmP8BfoMzG2UETtPoCdZaa4xZB/wkmKLPw/lLOTl4yjPAXcaYN3ES6Y+CxyJFSNdpjOkGvAY8a6296xSvUwjcYIy5AjiC839U24E1HVd6m4R0nUAlMAeoM8bcB/QHvoczyA9i5PNsYRLw5imOFxLZn2dIYuD+DEkM3J+hivb7s62i7v40xgwAXgLusdY+errzIune7HQtONbaJpxBetcAh3H6Qa+01h4yxswxxrRsJrsGpz+4HGftiVuttRuDzz0O/BNYhvMXcDMnDxZzVRuu8wacaYm3G2OOtfjz9+DzdwPv4fQFl+NMU51lrfWF83pOJ9TrDP6fxqU4U6wrgUU4n9+vgi8VK5/ncYXAqcYERPTneSaxdH+eSSzdn2cSS/fnmcTQ/XkHkAk80Orv4s8j9d6MCwQCH32WiIiISBTpdC04IiIiEvsUcERERCTmKOCIiIhIzFHAERERkZijgCMiIiIxRwFHREREYk5nXOhPRMLAGPMkzs7Bp3Mvzm7EbwOZ1tqwLLtvjIkHlgKft9ZuO8N5Hpw1SW6w1tpw1CYi7UctOCLSUb4O9Ar+mRE8NqHFsV/iLPbVC2dF03D5GrDuTOEGILjHzn18sLGliEQRLfQnIh3OGDMM2AD0t9budrGOFGAPcGGLlVU/6md24qzEurAjaxOR9qUuKhFxjTFmBi26qIwxAeC/gO/gbMi3Crge+BbOtgXVwHestc8Efz4TeAj4NBAA3gK+bq091TL4AJ8DjrQMN8aY7wNfAnoCW4DvWmv/0+Jn5uK0Ri1sh0sWkTBRF5WIRJqfAd/A2ZCwL85Gg9XAeODfwO+MMRnBc3+PE4QuBs7HCTmvGmNO9z9vlwOvHH9gjLkq+F7X4+xqvAD4pzGmS4ufeQX45BleU0QikAKOiESax6y1b1tr1+LsXnwMp1XF4my+mAr0N8YU4bTIXGetXRlslbkBZyPDS07z2uOATS0eFwKNQEmw6+w+4GrA2+KczUAGTgASkSih/yMRkUizo8X3dcBua+3xwYINwa/JQL/g99YY0/Ln03BadV46xWvnAhUtHv8FZ6ZXsTFmNTAP+LO1tr7FOZXBrzltvA4RcZFacEQk0nhbPfaf5ryE4LmjgVEt/gwC/nyan/EDcccfWGsPAWNxWnyWATcB64ODoo87/nvSF/IViIjrFHBEJFptARKBdGvtDmvtDqAUeBAn5JxKGc5gYgCMMVcDX7bWvmat/TpOy08NcFmLn+nZ4mdFJEqoi0pEopK11hpj5gFPG2PuAA4BP8EZnLz1ND+2GhjZ4nE88KAx5iDOjK1JQF7w++NGAlWc3HUmIhFOLTgiEs1uxAkjLwArgSxgprX2yGnOX4Az2woAa+0/gR/itPpsA34MfMVa+1aLn5kOvGKtVReVSBTRQn8i0mkYY9KA3cAl1to1IZzvAUpwZmot7uDyRKQdqQVHRDoNa20dTmvNHSH+yKeAYoUbkeijgCMinc2vgRGm1dzy1oKtN/cAt4WlKhFpV+qiEhERkZijFhwRERGJOQo4IiIiEnMUcERERCTmKOCIiIhIzFHAERERkZijgCMiIiIxRwFHREREYo4CjoiIiMScqN9N3BiTDIwHSgFthiciItJ5xAO9gJXW2saWT0R9wMEJN9onRkREpPOaBixpeSAWAk4pwF//+lfy8vLcrkVERETCpKysjDlz5kAwC7QUCwHHB5CXl0dBQYHbtYiIiEj4fWiIigYZi4iISMxRwBEREZGYo4AjIiIiMUcBR0RERGJOLAwyPiO/309FRQVHjhzB54uNZXJSUlIoKCggMTHR7VJEREQiUswHnH379hEXF0dhYSGJiYnExcW5XdJZCQQCVFZWsm/fPvr37+92OSIiIhEp5ruoamtryc/PJykpKerDDUBcXBzdu3enoaHB7VJERETOaNeBoyxdf4BAIBD29475FhwAjye2clwsBDUREYlNzT4/724oZcHSXWwqrgTg51+Zyrn9u4e1jk4RcERERKRjVdU08Np7Jfzn3d1UHnV6GVKTE7hkciGD+nYLez0KOCIiIvKxbdtTxfwlxSxZe4Bmnx+APrkZXD6liAvGFpCW4s6EGAUcERERaRNvs4/Faw/w0pJitu89AkBcHEwcmscVU/sz8pyerg+nUMARERGRkFQerec/y3bz6nslHDnWCEBGaiIXTezHZVP6k5ud5nKFH1DAiRDf/OY3yc7O5p577gHA5/Mxbdo0fv3rXzNx4kSXqxMRkc4qEAiweddhXlpSzLsbSvH5nRlRhb26cMXUIs4fk09KUuTFicirqIPd+8R7rNpyMGzvN25ILj/8wqSPPO+qq67i7rvv5u677yY+Pp6lS5eSkpLChAkTwlCliIjIyRq9Pha/v4/5i3dRfOAoAB5PHFNG9OaKqf0ZWtTd9W6oM+l0ASdSnXfeeXg8HpYvX855553H/PnzmTVrVkT/5RERkdhTXlV3ohuqpq4JgC7pSVw8qR+XTu5Pz26pLlcYmk4XcEJpTXGDx+Nh9uzZzJ8/n9GjR/PGG2/w73//2+2yRESkEwgEAmwsruSlJcW8t6GUYC8UAwuyuGJqEdNG5ZOUGO9ukW3U6QJOJLvqqqv43Oc+x6RJkxg0aJC2YhARkQ7V0NTMO2v289KSYnaXVgMQ74lj+sjezJpWhOnXLWp7EhRwIsiAAQPo168fDz/8MF/84hfdLkdERGJU+eE6Xl62i9eWl1BT5wWga0YyF0/ux6WTC+meFR3dUGeigBNhrrrqKh544AEuu+wyt0sREZEYEggE2LizkvlLilm+8YNuqHP6dGXWtCKmjuxNYkJ0dUOdiQJOhJkzZw5z5sxxuwwREYkRp+uGOn9UPrOm9cf0y3a5wo6hgCMiIhKDTtkNlZnMJZMKufS8QrK7pLhcYcdSwBEREYkRx2dDzV/84W6o2dOKmBJj3VBnooAjIiIS5Rq9Pt5Zs4/5iz/ohkqIj2PaiNjuhjoTBRwREZEo9cGifLs7ZTfUmXSKgOP3+/F4PG6X0W4CgYDbJYiIiEuO7w01f3Ex724sxR/sh+qM3VBnEtaAY4yZCfwMOAcoBx601v7OGJMM1ABNLU5fZq296GzfMz09nf3795Obm0tiYmLULlh0XCAQoLKykpSUzpvKRUQ6oyavj0Wt9oaK98Rx/uiCTtsNdSZhCzjGmD7Av4AbgReBscCrxpjdQCVw2Fqb197vW1BQQEVFBSUlJTQ3N7f3y7siJSWFgoICt8sQEZEwqDxaz4Klu3j1vRKqa512gKyMJC6ZXBgzi/J1hHC24BQCf7PWzg0+XmmMWQhMAfYDazviTT0eDzk5OeTk5HTEy4uIiLS7QCDA1t1VzF9SzLL1B/AFu6EGFGQxe1oRU0dG395Q4Ra2gGOtXQwsPv7YGJMNTAOeAS4Bcowx64FcYBHwDWvt/nDVJyIi4jZvs4/Faw8wf0kxO/YeAcDjiWPKyN7MnlbEkMLsqB9qES6uDDI2xmQB84DlON1V04ClwH2AF/g/YC4wwY36REREwqmqpoFXlu3m5Xd3c6SmEYDMtEQunlTIZef1p2c3dUO1VdgDjjFmEE6o2QzMsdb6gTtbnXMncMgY08dauzfcNYqIiITDjr1HmLd4J4vXHqDZ5wegsFcXZk0r4vwxBSSrG+pjC/csquk44eZx4LvW2kDw+H3A3621W4KnJgW/NoSzPhERkY7m8/l5d2Mp8xYVs2X3YQA8cTBpWB6zpw1g2IDu6oZqB+GcRTUAeAm4x1r7aKunRwDjjDHXBR8/Aiyw1h4KV30iIiIdqbq2iVff283LS3dRcdT5//f0lARmTuzH5VP6k9c93eUKY0s4W3DuADKBB4wxD7Q4/hhwK864mx3BmhYAXwpjbSIiIh2ipLSa+UuKeXv1Ppq8PgAKcjKYNa2IC8b2ITW5U6y5G3bhnEV1J63G2rQyJ1y1iIiIdCSfP8CqzWXMW1zM+h0VJ46PHZzD7GkDGDWoJx6PuqE6kmKjiIhIO6lr8PL6ij28tKSYsso6AFKS4vnE+L5cMbU/BTmZLlfYeSjgiIiInKUDh44xf0kxb67cQ32j0w2Vk53GrKn9+eSEfmSkJrpcYeejgCMiIvIxBAIB1m0/xLzFxazacpDj+yAPH9CDWdOKmDA0j3h1Q7lGAUdERKQNGpqaeWfNPuYtLmZPWQ0AiQkeZowpYNa0Ivr3znK5QgEFHBERkZBUHKnn5WW7eOXd3dTUeQHI7pLMZef155LJhWRlJLtboJxEAUdEROQMbMlh5i0qZmmLTS/P6dOV2dOKmDIyn8QEj8sVyqko4IiIiLTS7POzbP0B5i0qxu6pApxNL6eO7M3saQMYXNhNqw1HOAUcERGRoOOrDS9YuovK4GrDGamJXDypH5dPKdKml1FEAUdERDq9PWXVzFt88mrDfXIz/n97dx5e1XXee/wrCSQk5nkGjV4Y29hmMmYQkjzFI4NomsRJ45s8t2nqtEnTpLcZbhLH6c3kJM7QPE6aNk3dJmmNGAx4wDYSAszkCdtgXutIQmIehJg0oeHcP/aROVaEfAQ6OoN+n+fRY/Y++2y92/ss6dXa612L+xdlkT9zEgM023DM0R0TEZE+qa3Nz2t2gqdLy3n93UtLH7bPNnyzG63HUDFMCY6IiPQpjU0tvPTKQdZtKefwyToAUpKTKJg9mfsXZjJ5rGYbjgdKcEREpE84UVvPM9sqeW5HFXUNXpn3qGGp3LcggzvnTWVwWnKEI5SepARHRETilt/vx6pqWVtazstvHaUtUOZ9bfoI7l+UyfwbxpOUpDLveKQER0RE4k5Laxvb9hzh6S3lvFt9BoCkxAQW3zyJB3IzuWbK8AhHKOGmBEdEROLG+fqLPLf9/WXeg9P686Fb07lnfgajhqnMu69QgiMiIjHv4PHzrNtSwUuvHPzTMu9ZkxiQrF93fY3uuIiIxCS/388b755kbWk5r+4/8d7+mdPGsERl3n2eEhwREYkpTc2tlLx6iKe3lL+3mndy//Yy7wymjBsS4QglGijBERGRmHD6XCPPbKvk2e0HOFd3EYARQwZw7wJvNe8hA1XmLZcowRERkahWfugMT2+poPT1Q7S0emXe2ZOHsSQ3iwUzJmg1b+mUEhwREYk6rW1+du87xtrSct4urwEgMQFuvWE8S3KzmJ4xQuNrpEtKcEREJGo0NLXw4q5q1m2p4GiNt4xCako/7rxlKvctzGDcyIERjlBihRIcERGJuBO19azfWsnGHQeoa2wBYOyINO5flMkdc6eQNqB/hCOUWKMER0REImZ/1WnWbn7/MgrTM0awJDeLW64fT1KiHkPJlVGCIyIivaq1tY3tbx9l7eZy9lfVApeWUViyOJOcyVpGQa6eEhwREekVdQ3NvLCrinVbKjhR2wDAoNT+3DVvKvctzNQyCtKjlOCIiEhYHaupY92WCl7YVUVDDtoDTgAAIABJREFUk7eMwoRRA3kgN4vbZk9mQIp+FUnP06dKRER6nN/vZ1/ladaWlrPz7aMEhtcwI3sUSxZnMXvaWBI1vkbCqFcTHOfcHcD3gBzgBPBDM/uVcy4Z+AWwAmgFfmxm3+3N2ERE5Oq1tLaxbc8R1paWU3bwDAD9khLIu3kSS3KzyJw4NMIRSl/RawmOc24yUAR8ElgLzAKed84dAPIAB2QBQ4HnnHOHzew/eis+ERG5chcamtm44wDrtlRw6mwjAIPTkrlnfjr3LMhgxJABEY5Q+pre7MFJB35vZqsD27udcyXAAryk5yEzqwVqnXOPAZ8BlOCIiESxo6fqWLe1ghd2VtF40RtfM2nMIJbkZpE3axIDkjUSQiKj1z55ZrYF2NK+7ZwbASwCngTGA/uCDt8P3NBbsYmISOiCx9fsePso/sD4mptyRrNkcRYz3RiNr5GIi0hq7ZwbCjwN7AReDeyuDzqkHkjr7bhEROTyWlrbePnNI6zZ/P7xNYtneuNrMiZofI1Ej15PcJxz1+CNwdkHPAi0T3wQPAFCGnChl0MTEZFOeONrqli3tYJTZ7z5azS+RqJdb1dR5eIlN08AXzUzP9DonDuGN8j4cODQabz/kZWIiPSyYzV1PL2lgheD5q+ZOHoQSxdrfI1Ev96sosoC1gNfM7Ofd3j5SeCbzrk3gUHAl4Cf9lZsIiLi8fv97D9Qy+rNvvfNX3NjziiWLs7W+BqJGb2Zfj8MDAa+65wLnuPmn4FvAD8C9gKJwK/xenlERKQXtLa28fJb3vpQVu2tD9U+f83SxRpfI7EnpATHOXcDcDcwGxiDNxnfMWA3sN7MfB90DjP7IvDFLg55OPAlIiK9pL6xmY07q1m3pfx960PdPT+d+xZmanyNxKwuE5zAmJlHgPnALrxxMT4gCRgFfBz4gXNuM/ComZWGN1wREekJJ2rrWbelgo07q6hvbAG89aGWLM6iYJbWh5LYd9lPsHPu34Dr8JZQWGZmZy5z3BDgo8Djzrk3zeyhcAQqIiJXr+xgLWtKytn65hHaAgNsrs8aydLcLOZMH6fxNRI3ukrRN5jZpz7oBGZ2DvgV8Cvn3Ioei0xERHpEa5uf3fuOsWZzOXsragBITExg8c2TWLI4k5zJwyMcoUjPu2yCY2ZF3T2Zma28unBERKSnNF5s4aXdB3m6tJwjp+oASBvQj7vmpXP/wkxGD0/9gDOIxK5QBxmnAf8H+E8zK3PO/Qpvkr5dwINmdjSMMYqISDfUnmtk/bZKnn25kvP1zQCMGZ7KA7lZ3DF3CmkD+kc4QpHwC3UU2U+BAuB/nHNL8BbH/BtgKfAz4M/CE56IiISq6ug51mwup+S1Q7S0tgFwzZRhLF2czfwbxpOUlBjhCEV6T6gJzhLgPjPb65z7CvCCmf2Lc24bsD184YmISFf8fj9vvHuSNZvLec1OAJCQAPOuH8fSxdlMzxhBQoIGDkvfE2qCkwocd84lAncB3wrs9+PNiSMiIr2ouaWN0tcPsWZzOQeOngMgJTmJ2+dM4YFFmUwYPSjCEYpEVqgJzm68MTgngeHAaufcBOBRYEeYYhMRkQ4u1F/k2e0HWL+1gtPnmgAYPjiF+xZmcvf8dAanJUc2QJEoEWqC8zng90A68NdmdsQ593O8BTKXhSk2EREJOFZTx9rScl7cVU3jRa/jfOq4wSxdnM3imRPp3y8pwhGKRJeuJvqbD+wwszYz2wfc1OGQr5rZ+bBGJyLSx+2vOs3qEh873rq08OVN14xm2eJsbnajNb5G5DK66sF5DHDOuS3ARmBj8JpTSm5ERMKjtc3Prr1HWV1SzjsHTgNa+FKku7qa6G9+YBmG24A7gC8655KAF/ASnpfMrLZ3whQRiX/tE/Ot3VzO0RpvYr6Bqf25+9Z07luYwcihmphPJFRdjsEJLMOwOvCFcy4DuBP4CPCEc64cr2T86+EOVEQkXtWeb2TD1kqeCZqYb+yINB7IzeSOuVNJ1cKXIt3WrVZjZpVcWncqEZiD17sjIiLdVH3s0sR8zS2XJuZblpfNrddrYj6RqxFyguOcy8NbXTylw0v1PRmQiEg88/v9vF1ew6oSH6+8cxzwJua75bpxLMvTxHwiPSXUtagex1uaoRpo7PCyH/hxD8clIhJXWlvb2PbmEVaX+PAdOgtAcr9EbpszhSWLs5ioiflEelSoPTh/AXzKzH4XzmBEROJNfWMzL+yqZm1pOSdrGwAYOiiZe+dncM+CDIYO6tgpLiI9IdQEpx5v5XAREQlBzdkG1m2p4LntB6hrbAFg4uiBLF2cTf7syaT018R8IuEUaoLzHeAx59znAgONRUSkEweOnmN1iY/S1w/R0urNzHdd5kiWLc5izvRxJCZqfI1Ibwg1wXkH+H+Azzn3Jy+amf4UEZE+y+/3s6fsJKtLLq3onZgAC26cwPK8bK6ZMjzCEYr0PaEmOL/GW1Tzt6hqSkQEgJbWNra+cZjVJeVUHPEGDqckJ3HH3Cksyc1i3MiBEY5QpO8KNcGZDNxtZhXhDEZEJBbUNzazcWcVa0srOHXGGzg8bHAK92tFb5GoEWqC8wKQCyjBEZE+q33g8LPbD1AfGDg8acwgluVlkzdzEskaOCwSNUJNcHYCv3DOFQI+oDn4RTP7h54OTEQkWnQ2cPj6rJEsy8tm9rSxGjgsEoVCTXDuAHYDg4CbOrzm79GIRESigN/v503fKVaV+HhtvwYOi8SakBIcM8sPdyAiItGgtbWNrXuOsHqzj/JDQQOHAzMOa+CwSGy4bILjnPsm8AMzawjlRM65wcCXzewbPRWciEhvaWhq4YWdVawtLedE0IzD9y3M5J75GQwZqIHDIrGkqx6cs8Be59xKYJWZ7eh4gHMuAZgNfBxYDvwklG/qnJsLrDezMYHtFOA8cDHosJfN7M6QrkJE5ArVnmtk3dYKnnn5AHUN3vBCzTgsEvsum+CY2eOB5OYfgI3OuRa8Cf9OAQnAaLzVxROAfwcWmFl1V98skBB9Gnisw0s3AKfNbNwVXoeISLccPH6eNZvL2fTKQVpa2wC4Nn0Ey/KyueU6zTgsEuu6HINjZoeAv3XOfQXIA2YBY4E2vMqqR4BiM2sK8fs9AtyLt/TD14P2zwLe6FbkIiLd5Pf72Vd5mtUlPnbuPQZAQgLMu34cy/NyuDZjRIQjFJGeEuog4zpgQ+DrajxhZt9wzuV12D8TGOOcexMvgSoFvmBmh6/y+4mI0NrmZ+fbR1lV4sOqagHo3y+R2+ZMYeniLCaOHhThCEWkp4VaJt4jzOzIZV6qA7YB38abY+dnwGpgbi+FJiJxqKm5lU2vHGRNiY8jp+oAGJzWn3sWZHDfgkyGDU6JcIQiEi69muBcjpl9MXjbOfdF4KRzbrKZHYxQWCISo87XX+SZbZWs21rB2Qte7cKYEWkszc3ijrlTGJASFT/6RCSMoqKVO+e+DfzBzN4J7Gqvx2yMUEgiEoOOn65nbWk5G3dW0XSxFYDsSUNZnpfD/BnjSUpKjHCEItJbQkpwnHPpZnYgjHHMAGY75z4W2P4psMHMTobxe4pInCg/dIZVJT627jlCW5s3ufpMN4bl+dnMyB5FQoIqokT6mlB7cHzOuR3AfwJPmVlND8fxabxxN75ATBuAv+zh7yEiccTv9/P6uydZXezjjTLvb6GkxATyZk1ieV42GROGRjhCEYmkUBOcTOBjwGeBnzrnNgL/BawNdabjYGZWAgwL2q4BHuzueUSk72kJLKWwqriMyiPnAEhNSeLOW9J5IDeTMcPTIhyhiESDUMvEq4HvAd9zzl0PfAT4CvBr59xq4EkzezF8YYpIX9fZUgrDBqfwwKJM7r41nUFpWkpBRC65kkHGh4ByoALIwVtd/A7n3HngITPb3oPxiUgfd+Z8E+u3VrBhWyUXgpZSWJaXQ/6sSSRrKQUR6USog4wHAkvwem7uBI4Dvwe+bmZ7nXOJwC+B/wamhClWEelDjpy6wJrN5by0q5qLLd5SCtOmDqewIIe507WUgoh0LdQenBN4C2GuAj4UGEPzHjNrC4zLWdSz4YlIX/NudS2rSnxsf/MIgYIobrluHMvzs5meMTKywYlIzAg1wXkIeLqzNaecc2PM7ISZrcJLgEREusXv9/OanWBVsY83facA6JeUwG2zJrMsL5vJYwdHOEIRiTWhJjh/BMYB75uXxjk3BdgHaCEXEem2ltY2tr5xmKJiHweOehVRaQP6cfet6dy/KJORQ1MjHKGIxKrLJjjOuY8CywKbCcBvnHMde3CmAqfDFJuIxKnGphY27qpizeZyTgYqokYMSeGBRVl86NZ0Bqb2j3CEIhLruurBeQG4Ay+5AWgIfLXzAzuBfw9LZCISd85eaGLd1gqe2VbJ+XqvImrSmEEsz8smb9Yk+vdTRZSI9IzLJjhmdgr4FIBz7gDwQzOr752wRCSeHKupY3WJjxdVESUivaSrR1T3AC+YWTOwG8hzznV6rJk9E57wRCSW+Q6dYVWxj217Dr9XETVn+lgK83O4LlMVUSISPl09olqPN7D4RODfl+MH1K8sIoBXEbWn7CRFm96/RlTBrEksz89m6rghEY5QRPqCrh5RJXb2bxGRzrS2+Xn5TW+NKN+hs8ClNaKW5GYxergqokSk94S8VINz7lPAOTNbGdj+H2CdmT0ZruBEJPo1NbeyaXc1q0vKOVpTB8CwQSncvyiTe+ZrjSgRiYxQl2r4GvD3eKuJt3sLeNw5N9LMHg9HcCISvS7UX+SZlw+wbksFZy54M0iMHzmQZXlZFMyZQorWiBKRCAq1B+czwEfMbGP7DjN71Dn3Ct4aVEpwRPqImrMNrNlczvM7DtDQ1ApA1qShFObnMH/GBJJUESUiUSDUBGc4UNXJ/nJgbM+FIyLR6uDx86wq9lHy2kFaWr2SqJtyRrOiIIcZOaNISFBiIyLRI9QEZwfwj865/21mLQDOuSS8x1a7wxWciETe/qrTFG0qY+feY/j9kJgAC2+cQGF+DtmTh0U6PBGRToWa4HwJeAmods69iVcafkPg/XeHKTYRiRC/38+r+0+wclMZeytqAOjfL5Hb5kxhWV4WE0Zp+TkRiW4hJThmtsd5s/x9BLgWuAisBf7LzM6HMT4R6UWtrW1s6bD45cAB/bhnQQb3L8xk+JABEY5QRCQ0IZeJm1mNc+554CCQCOxXciMSHxovtvDirmpWby7nxGlvRZYRQ1JYkustfpk2QItfikhsCbVMfBDwr8AKoBlvAc5+zrkXgEIzqwtfiCISLufqLrJhWyXrt1Zwru4iABNHD2R5fg75WvxSRGJYqD04P8Ybc3MrlwYVzwF+A3wf+FzPhyYi4XKytoE1pT427qii8aJX6n3NlGEU5udwy/XjVeotIjEv1ARnObDMzHYF7dvlnHsYWIkSHJGYUH3sHEXFPja/dojWwOqXM90YVhTkcH3WSJV6i0jcCDXBSQROdbL/NKByCpEot6+yhqJNPnbtOwZ4pd65N0+kMD+HzIlDIxydiEjPCzXBKQW+5Zz7hJldBHDOpQDfBLaEKzgRuXJtbX5e2X+cok1l7Ks8DUByv0RunzuFZXnZjBs5MMIRioiET3fmwdkKHHTOvRHYdyPQCHwoHIGJyJVpaW2j9PXDFBWXUX3MK3QcmNqf+xZkcN/CTIYNTolwhCIi4RfqPDg+59y1wIPAdKABeApvHpyGMMYnIiFqbGph464q1mwu52St1yxHDh3A0sVZ3HnLVJV6i0if0p15cGqBX4QxFhG5AmcvNAVKvSs5X++Vek8aM4jC/GwWz5xM/36JEY5QRKT3XTbBcc7txluS4QOZ2dzufFPn3FxgvZmNCWwn4yVPK4BW4Mdm9t3unFOkrzlxup41peVs3FlFU6DU200dzoqCHOZOH0eiSr1FpA/rqgdnfU9/M+dcAvBp4LEOLz0COCALGAo855w7bGb/0dMxiMS6qqPnKCouo/T1w++Ves+a5pV6X5epUm8REegiwTGzR8Lw/R4B7gW+A3w9aP8ngYcCj8FqnXOPAZ8BlOCIBOytqKGouIzd+44DkJiYwOKbJ1FYkE3GBJV6i4gEC3kMjnPuw8CXgRxgJvDXwDEz69gb05UnzOwbzrm8oPMOA8YD+4KO2483c7JIn9bW5mf3vmMUFft450Cg1Lt/EnfOncLSvGzGjkiLcIQiItEp1LWoHsJ7rPQj4P8Gdu8Hfuyc62dm3wvlPGZ2pJPd7RMF1gftqwf0k1v6rOaWNkpfP0RRsY+Dx71S70Gp/bl3obeq99BBKvUWEelKqD04fw981syecs59FcDMfuOcq8VLfEJKcC6jfaHO1KB9acCFqzinSExqaGrh+R1VrN3s49TZRgBGDR3AksXZ3DVvKqkpIXe6ioj0aaH+tMwCXulk/xvAuKsJwMxqnXPH8AYZHw7snsb7H1mJxLWzF5pYt7WCDVsrudDQDMDksYMpzM8m9+ZJKvUWEemmUBMcA24H/qXD/g/jPaq6Wk8C33TOvYn3yOpLwE974LwiUe346XrWlPjYuKuai81eqfe16SMozM9mjkq9RUSuWKgJzleBlc652YH3/JVzLhu4D2/umqv1DbzxPXvxFvb8NfBED5xXJCpVHjlL0SYfW/Ycpi1Q6j1n+lgK871SbxERuTqhLtXwbGByvi8DbwN3AO8A88zste5+UzMrAYYFbTcCDwe+ROKS3+/n7YoaijaV8er+E4BX6p0/axLL83NIHz8kwhGKiMSPrmYyvgd4zszaAMxsL/BQL8UlEjfa2vzs3HuMouIyrKoWgJTkJO68ZSpLc7MYo1JvEZEe11UPzlqgxjn3B+B3ZvZGF8eKSAfNLW2UvHqQVSU+Dp3wigIHp/XnvoWZ3LsgQ6XeIiJh1FWCMxH4SODr8865vcDvgN9fZj4bEQHqG5u9Uu/ScmraS72HpbIssKr3AJV6i4iEXVdLNZwAfgb8zDmXDnwU+ATwXedcMd4yCkVm1tAbgYpEuzPnA6Xe2yqpC5R6Tx03mOX5OeTePJF+SSr1FhHpLaEOMj4AfBcvuZkOfAyvsuqXzrkiM/tf4QtRJLodq6ljVYmPl3ZVc7GlDYDpGSNYUZDD7GvHavFLEZEI6HZfuZntCyyGacAX8Hp1lOBIn1Nx+CxFm8rYuucwgUpv5k4fx4qCHK7NGBHZ4ERE+rjuLLY5FFiKN7nfbUAF8F/A8vCEJhJ9/H4/b5WfYuVLZbz+7kkAkhITKJg1ieX52Uwdp1JvEZFo0GWCE5TU/BneTMbngD8C3zKz3eEPTyQ6tLX52fH2UYqKy3i3+gwAA5KTuHPeVJbmZjN6eOoHnEFERHpTV/PgbMDrqWkD1uHNWPycmbX0UmwiEdfc0krxq4dYVVzG4ZPeurBDBiZz/yKv1HtwWnKEIxQRkc501YMzCG9m4afM7FwvxSMSFeobm3luu1fqffqcV+o9Zngqy/KyuX3uFAYkq9RbRCSadVUmvrg3AxGJBrXnG1m3pYJntlVS1+h1VqaPH0JhfjYLb1Kpt4hIrNCfoSLA0VN1rC7x8eLuapoDpd7XZY5kRUEOs6aNUam3iEiMUYIjfVr5oTMUFfvYFlTqfct1Xqn3tHSVeouIxColONLn+P1+3vSdYuWmMt4IlHr3S0qgYOZkludnM3ns4AhHKCIiV0sJjvQZrW1+drx1lJXFZfgOeqXeqSlJ3DUvnSW5WYwaplJvEZF4oQRH4t7F5laKXz3IqmIfR055pd5DByVz/8JM7lGpt4hIXFKCI3GrrqGZZ7cf4OnScmrPNwEwZkQaywOl3in9kyIboIiIhI0SHIk7p8818nRpOc9uP0B9oNQ7Y8IQCvNzWHjjBJJU6i0iEveU4EjcOHLygreq9+6DtLR6pd43ZI1iRUEON7vRKvUWEelDlOBIzCs7WEvRJh8vv3UEvx8SEuDWG8azoiCHa6YMj3R4IiISAUpwJCb5/X7eePckRcVl7Ck7BXil3vmzvFLvSWNU6i0i0pcpwZGY0trm5+U9R1hZXEbF4bMApKb04+5b03kgN5ORQ1XqLSIiSnAkRlxsbuWl3dWsLinnaI1X6j1scAoPLMrk7vkZDErtH+EIRUQkmijBkah2oaGZZ1+u5OktFZwJlHqPHzmQZfnZ3DZ7Mskq9RYRkU4owZGoVHO2gbWlFTy3/QANTV6pd9akoRTm5zB/xgSSElURJSIil6cER6LKoRPnWVXso/jVQ++Vet+YM4rC/Bxuukal3iIiEholOBIV3q2uZeWmMna8ffS9Uu8FN06gMD+bnMkq9RYRke5RgiMR4/f7ed1OsnJTGW+Vt5d6J3LbnMksz8tmwuhBEY5QRERiVdQkOM65TwG/ApqCdj9sZr+LUEgSJq2tbWx78whFm3xUHPFKvdMGtJd6ZzFiyIAIRygiIrEuahIcYCbwIzP7x0gHIuHRFCj1XlXs4/jpegCGD07hgdws7r41nYEq9RYRkR4STQnOLOCnkQ5Cet6F+otseLmSdVsqOHvhIgDjRw1keV42BSr1FhGRMIiKBMc5lwTMAD7hnPsxUA/8Bvi+mfkjGpxcsZqzDazZXM7zOw7Q0NQKeKXeKwpyuPUGlXqLiEj4REWCA4wGXgF+BywHrgXWAueAX0YwLrkCB497pd4lrx2kpdXLT2/KGc2Kghxm5IxSqbeIiIRdVCQ4ZnYMWBy06w3n3M+BQpTgxIz9Vacp2lTGzr3H8PshMVDqvSI/h+zJwyIdnoiI9CFRkeA4564DPmxm3wzanQw0RigkCZHf7+fV/ScoKi7j7fIaAPr3S+S2OVNYlpfFhFEq9RYRkd4XFQkOcAb4e+fcIeBfgZuBvwU+F9Go5LJaW9vYsucIRZvKOHD0HOCVet8zP4MHFmUyXKXeIiISQVGR4JjZYefcA8APgJ8Ap4BHzWxlZCOTjhovtvDirmpWby7nRKDUe8SQFJbkZvGhW9NJG6BSbxERibyoSHAAzGwTMDvScUjnztdfZMM2r9T7XJ1X6j1x9ECW5eVQMHsS/fup1FtERKJH1CQ4Ep1O1jawptTHxh1VNF70Sr1zJg+jsCCHedePV6m3iIhEJSU40qnqY+coKvax+bVDtLZ5pd4z3RgKC7K5IUul3iIiEt2U4Mj7vFN5mqJir9QbvFLv3Jsmsjw/m6xJKvUWEZHYoARH8Pv9vPLOcVZuKmNf5WkAkvslctvcKSzPy2bcyIERjlBERKR7lOD0YS2tbZS+fphVxWVUHTsPwMDU/ty7IIP7F2YybHBKhCMUERG5Mkpw+qDGphY27qpizeZyTtY2ADBy6ACW5GZx17ypKvUWEZGYpwSnDzl7oYkN2ypZv7WS8/Xtpd6DWFGQzeKZk+nfLzHCEYqIiPQMJTh9wInT9awpLWfjziqaAqXebupwCvNzuOW6cSSq1FtEROKMEpw4VnX0HEXFZWx+/TBtgVLvWdPGsKIgh+syR6rUW0RE4pYSnDi0t6KGouIydu87DkBiYgJ5MyexPD+bjAlDIxydiIhI+CnBiRNtbZdKvd85ECj17p/EnXOnsDQvm7Ej0iIcoYiISO9RghPjvFLvQxQV+6gOlHoPSu3PvQu9Uu+hg1TqLSIifY8SnBjV0NTCxp1eqfepM16p96ihA1iyOJu75k0lNUW3VkRE+i79FowxZy80sX5rJRu2VXC+vhmAyWMHsTwvh8UzJ6nUW0REBCU4MeP46XrWbPaxcWc1F5u9Uu9pU4ezoiCHOdNV6i0iIhJMCU6UqzxyllXFPkrfuFTqPfvasawoyGF6xgiVeouIiHRCCU4U8vv97K2oYeWmMl7dfwIIlHrPmkRhfg7p44dEOEIREZHopgQnirS1+dm17xgrN5VhVbWAV+p917ypLM3NYoxKvUVEREKiBCcKNLe0sfm1gxQV+zh04gIAg9P6c9/CTO5dkKFSbxERkW5SghNB9Y3N75V615xtBGDUsFSW5WVx59ypDFCpt4iIyBXRb9AIOHO+iXVbK9iwrZK6Bq/Ue8q4wRTm55B780T6JanUW0RE5GoowelFx2rqWF3i48Vd1VxsaQNgesYICgtymD1trEq9RUREeogSnF5QcfgsRZvK2LrnMIFKb+ZOH0dhQTbTM0ZGNjgREZE4pAQnTPx+P2+Vn6Jok4/XzCv1TkpMoGCWt6r31HEq9RYREQkXJTg9rLXNz863j1JUXMa71WcAGJCcxJ3zprIkN4sxw1XqLSIiEm5KcHpIc0srm145xOqSMg6frANgcFoy9y/ySr2HDEyOcIQiIiJ9hxKcq1Tf2Mxz2w+wtrSc0+eaABgzPJWli7O5Y+4UlXqLiIhEgH77XqHa842s21LBM9sqqWtsASB9/BAK87NZeJNKvUVERCIpahIc59yNwBPADKAC+JSZ7Y5sVH/q6KlAqffuapoDpd7XZY5kRUEOs6aN0eKXIiIiUSAqEhznXDKwFngcyAUKgY3Oualmdi6iwQWUHzpDUbGPbUGl3rdcN47C/ByuzRgR2eBERETkfaIiwQHygP5m9nhg+4/Ouc8Bfw78S8SiAvZW1PDfLxivv3sS8Eq9b5s9ieV52UxRqbeIiEhUipYEZzrwTod9+4EbIhDLe46cvMBXf7mVNr9X6n3XvHSW5GYxenhqJMMSERGRDxAtCc4goL7DvnogopPGDB8ygNvnTmXMiFTumZ/B4DSVeouIiMSCaElw6oCO3SJpwIUIxPKe1JR+/M2Hb4pkCCIiInIFoqWWeR/gOuybFtgvIiIi0i3R0oNTDCQ45/4O+AVeFdUMYHVEoxIREZGYFBU9OGZ2EbgbL7E5DXwNWGpmJyMamIiIiMSkaOnBwczeBhZGOg4RERGJfVHRgyMiIiLSk5TgiIiISNxRgiMiIiJxJ2rG4FyFJIBjx45FOg4RERHpRUG/+5M6vhYPCc54gAcffDDScYiIiEhkjAfKg3fEQ4KzG1gEHAVaIxyLiIiI9J6oLxbZAAAJDUlEQVQkvORmd8cXEvx+f++HIyIiIhJGGmQsIiIicUcJjoiIiMQdJTgiIiISd5TgiIiISNxRgiMiIiJxRwmOiIiIxB0lOCIiIhJ3lOCIiIhI3ImHmYy7zTl3I/AEMAOoAD5lZn8yC6Jzbgrwr8A84ATwN2b2TOC1BOBR4C+BZOC3wJfNrKVXLiIE3bjOWcBPAsedA34DPGpm/sDr1cBIoH1WyMNm5sJ/BaHpxnUWAC8ADUG7v29mj8bL/XTOLQKe7fDWFKDSzK4JHBPV97Odc24usN7Mxlzm9Zhun+1CuM6Ybp/tQrjOmG6f7bq6zlhvn865O4DvATl4be6HZvarTo6LirbZ5xIc51wysBZ4HMgFCoGNzrmpZnauw+F/BLYD9wILgTXOuZvMrALv5iwHZgJNwGrgq8C3e+VCPkCo1+mcSwM2AN8B8oFM4HngGPBr59woYCIwxMzqevcqPlg37+dM4Ckz+0gnp4qL+2lmW4BBQe+bjDeF+ecC21F9P+G9H4CfBh77gENjtn1CaNcZ6+0TunU/Y7Z9QmjXGcvtMxBrEfBJvJ9Fs4DnnXMHzOz5DodHRdvsi4+o8oD+Zva4mTWb2R+BvcCfBx/knLsGmA18w8wumtkm4Gm8DzB4N/lxMztkZieBbwGf6aVrCEUeIVwnMBnYbma/MLNWMysD1uB9KMH7EJdFW2MLkkdo1wnetbxxmfPEy/3s6N+AJ81sY2A72u8nwCPAZ/F+qXcqDtonhHCdxH77hNCuE2K7fULo1xksltpnOvB7M1ttZm2B3uMSYEHwQdHUNvtigjMdeKfDvv3ADZ0cV93hgxZ83HRgX4fXJjjnRvRgrFcjpOs0z7L27UBPwd3A64FdM4FE59wu59xJ59zzzrlrwxh3d4V6P8G7ltudc1XOuWrn3A+dcylB54n5+xnMObc08L5vBO2O9vsJ8ISZzQJe6eKYWG+fEMJ1xkH7hNDuJ8R2+4TQrxOIvfZpZlvM7K/atwP/7xdx6bPYLmraZl9McAYB9R321QNp3Tyu4+vt/+54nkgJ9TrfE/hh8ofAcU8EdrcCu/C6FKfifZifDXSdR4OQrtM51w84hNcdei1QANyO9yy4s/PE/P0EvgZ8z8yCxzRE+/3EzI6EcFist89Qr/M9Mdo+Q7rOOGif3b6fxGj7BHDODcXrldmJ97gqWNS0zT43BgeoA1I77EsDLnTzuI6vt9+cjueJlFCvEwDn3Di856ttwO3tjc7MftDhuK8Af43Xlbqlh2O+EiFdZ2AA221Bu3zOuX8Cvg/8QyfnifX7OQO4Dvhd8P4YuJ+hivX22S0x3D5DEgfts1tiuX0GHkGtxeuFedDM2jocEjVtsy/24OwDOo5In8b7u8zaj5vinEu9zHEdzzMNOGpmZ3ow1qsR6nXinJuON9DNh/fDszbotS845xYGHZ6Elxg39njEVyak63TOTXTOPRbo4m+XzKXriJv7GbAEeLbjQOsYuJ+hivX2GbIYb58hiYP22V0x2T6dc7l4vTZrgBVm1llcUdM2+2IPTjGQ4Jz7O+AXeNUoM/C6Rt9jZuac2wP8UyCLno/3obw1cMiTwJeccy/hZaTfCuyLFiFdp3NuOLAR+KOZfamT86QDn3DO3QecwfuLqgx4LXyhd0tI1wnUAA8C9c65bwMZwNfxBvlBnNzPIPOAlzrZn05038+QxEH7DEkctM9QxXr77K6Ya5/OuSxgPfA1M/v55Y6LprbZ53pwzOwi3iC9QuA03nPQpWZ20jn3oHMuuJusEO958Am8uSc+bWZvB157AngKeBnvA7iP9w8Wi6huXOcn8MoSP+ucuxD09YfA6/8I7MB7FnwCr0z1fjNr7c3ruZxQrzPwl8bdeCXWNUAp3v37ceBU8XI/26UDnY0JiOr72ZV4ap9diaf22ZV4ap9diaP2+TAwGPhuh8/i96O1bSb4/f4PPkpEREQkhvS5HhwRERGJf0pwREREJO4owREREZG4owRHRERE4o4SHBEREYk7SnBEREQk7vTFif5EpBc45/4db+Xgy3kEbzXiYmCwmfXKtPvOuSRgG/AXZvZuF8cl4s1J8gkzs96ITUR6jnpwRCRcPg+MD3zlBfbNDdr3GN5kX+PxZjTtLX8L7OkquQEIrLHzbS4tbCkiMUQT/YlI2DnnrgfeAjLM7EAE4xgAVAMFQTOrftB7yvFmYi0JZ2wi0rP0iEpEIsY5l0fQIyrnnB/4KPAVvAX5XgE+DnwZb9mCc8BXzOzJwPsHAz8CVgB+YBPweTPrbBp8gI8AZ4KTG+fc/wX+EhgNvAN81cyeDXrParzeqJIeuGQR6SV6RCUi0eZ7wBfwFiScgrfQ4DlgDrAK+JVzblDg2F/jJUJ3AYvxkpznnXOX++PtXuC59g3n3LLA9/o43qrGG4CnnHNDgt7zHHB7F+cUkSikBEdEos0/m1mxmb2Bt3rxBbxeFcNbfDEVyHDOZeL1yHzMzHYHemU+gbeQ4Ycuc+7ZwN6g7XSgCagKPDr7NrAcaA46Zh8wCC8BEpEYob9IRCTa+IL+XQ8cMLP2wYKNgf+mAFMD/zbnXPD70/B6ddZ3cu6xwKmg7f/Eq/SqcM69CjwN/NbMGoKOqQn8d0w3r0NEIkg9OCISbZo7bLdd5rh+gWNvBm4K+roG+O1l3tMGJLRvmNlJYBZej8/LwEPAm4FB0e3af062hnwFIhJxSnBEJFa9A/QHBpqZz8x8wFHgh3hJTmeO4Q0mBsA5txz4jJltNLPP4/X8nAfuCXrP6KD3ikiM0CMqEYlJZmbOuaeB/3DOPQycBP4Jb3Dy/su87VXgxqDtJOCHzrnjeBVb84BxgX+3uxGo5f2PzkQkyqkHR0Ri2SfxkpE1wG5gKHCHmZ25zPEb8KqtADCzp4Bv4vX6vAt8B/icmW0Kek8u8JyZ6RGVSAzRRH8i0mc459KAA8CHzOy1EI5PBKrwKrW2hDk8EelB6sERkT7DzOrxemseDvEtS4AKJTcisUcJjoj0NT8BZrgOteUdBXpvvgb8Va9EJSI9So+oREREJO6oB0dERETijhIcERERiTtKcERERCTuKMERERGRuKMER0REROLO/wdkT47KCCgaNwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 576x576 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "subplot(2, 1, 1)\n",
    "plot_position(results)\n",
    "\n",
    "subplot(2, 1, 2)\n",
    "plot_velocity(results)\n",
    "\n",
    "savefig('figs/chap11-fig02.pdf')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "According to this model, we should be able to make this run in just over 2 seconds."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "2.010095574667565 second"
      ],
      "text/latex": [
       "$2.010095574667565\\ \\mathrm{second}$"
      ],
      "text/plain": [
       "2.010095574667565 <Unit('second')>"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "t_final = get_last_label(results) * s"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "At the end of the run, the car has gone about 28 meters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>values</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>x</th>\n",
       "      <td>28.385825236293645 meter</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>v</th>\n",
       "      <td>27.777777777777775 meter / second</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "x             28.385825236293645 meter\n",
       "v    27.777777777777775 meter / second\n",
       "Name: 2.010095574667565, dtype: object"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "state = results.last_row()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we send the final state back to the slope function, we can see that the final acceleration is about 13 $m/s^2$, which is about 1.3 times the acceleration of gravity."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "99.99999999999999 kilometer/hour"
      ],
      "text/latex": [
       "$99.99999999999999\\ \\frac{\\mathrm{kilometer}}{\\mathrm{hour}}$"
      ],
      "text/plain": [
       "99.99999999999999 <Unit('kilometer / hour')>"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "v, a = slope_func(state, 0, system)\n",
    "v.to(km/hour)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "12.918946212080874 newton/kilogram"
      ],
      "text/latex": [
       "$12.918946212080874\\ \\frac{\\mathrm{newton}}{\\mathrm{kilogram}}$"
      ],
      "text/plain": [
       "12.918946212080874 <Unit('newton / kilogram')>"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "a"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "1.3182598175592728 dimensionless"
      ],
      "text/latex": [
       "$1.3182598175592728\\ dimensionless$"
      ],
      "text/plain": [
       "1.3182598175592728 <Unit('dimensionless')>"
      ]
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "g = 9.8 * m/s**2\n",
    "(a / g).to(UNITS.dimensionless)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It's not easy for a vehicle to accelerate faster than `g`, because that implies a coefficient of friction between the wheels and the road surface that's greater than 1.  But racing tires on dry asphalt can do that; the OEM team at Olin has tested their tires and found a peak coefficient near 1.5.\n",
    "\n",
    "So it's possible that our no slip assumption is valid, but only under ideal conditions, where weight is distributed equally on four tires, and all tires are driving."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise:** How much time do we lose because maximum torque decreases as motor speed increases?  Run the model again with no drop off in torque and see how much time it saves."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Drag"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this section we'll see how much effect drag has on the results.\n",
    "\n",
    "Here's a function to compute drag force, as we saw in Chapter 21."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [],
   "source": [
    "def drag_force(v, system):\n",
    "    \"\"\"Computes drag force in the opposite direction of `v`.\n",
    "    \n",
    "    v: velocity\n",
    "    system: System object\n",
    "    \n",
    "    returns: drag force\n",
    "    \"\"\"\n",
    "    rho, C_d, area = system.rho, system.C_d, system.area\n",
    "    \n",
    "    f_drag = -np.sign(v) * rho * v**2 * C_d * area / 2\n",
    "    return f_drag"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can test it with a velocity of 20 m/s."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "-72.0 kilogram meter/second<sup>2</sup>"
      ],
      "text/latex": [
       "$-72.0\\ \\frac{\\mathrm{kilogram} \\cdot \\mathrm{meter}}{\\mathrm{second}^{2}}$"
      ],
      "text/plain": [
       "-72.0 <Unit('kilogram * meter / second ** 2')>"
      ]
     },
     "execution_count": 34,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "drag_force(20 * m/s, system)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's the resulting acceleration of the vehicle due to drag.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "-0.24 meter/second<sup>2</sup>"
      ],
      "text/latex": [
       "$-0.24\\ \\frac{\\mathrm{meter}}{\\mathrm{second}^{2}}$"
      ],
      "text/plain": [
       "-0.24 <Unit('meter / second ** 2')>"
      ]
     },
     "execution_count": 35,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "drag_force(20 * m/s, system) / system.mass"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can see that the effect of drag is not huge, compared to the acceleration we computed in the previous section, but it is not negligible.\n",
    "\n",
    "Here's a modified slope function that takes drag into account."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [],
   "source": [
    "def slope_func2(state, t, system):\n",
    "    \"\"\"Computes the derivatives of the state variables.\n",
    "    \n",
    "    state: State object\n",
    "    t: time\n",
    "    system: System object \n",
    "    \n",
    "    returns: sequence of derivatives\n",
    "    \"\"\"\n",
    "    x, v = state\n",
    "    r_wheel, speed_ratio = system.r_wheel, system.speed_ratio\n",
    "    mass = system.mass\n",
    "    \n",
    "    omega2 = v / r_wheel * radian\n",
    "    omega1 = omega2 / speed_ratio\n",
    "    tau1 = compute_torque(omega1, system)\n",
    "    tau2 = tau1 / speed_ratio\n",
    "    F = tau2 / r_wheel\n",
    "    a_motor = F / mass\n",
    "    a_drag = drag_force(v, system) / mass\n",
    "    \n",
    "    a = a_motor + a_drag\n",
    "    return v, a "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And here's the next run."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>values</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>success</th>\n",
       "      <td>True</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>message</th>\n",
       "      <td>A termination event occurred.</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "success                             True\n",
       "message    A termination event occurred.\n",
       "dtype: object"
      ]
     },
     "execution_count": 37,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "results2, details = run_ode_solver(system, slope_func2, events=event_func)\n",
    "details"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The time to reach 100 kph is a bit higher."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "2.034354320417196 second"
      ],
      "text/latex": [
       "$2.034354320417196\\ \\mathrm{second}$"
      ],
      "text/plain": [
       "2.034354320417196 <Unit('second')>"
      ]
     },
     "execution_count": 38,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "t_final2 = get_last_label(results2) * s"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "But the total effect of drag is only about 2/100 seconds."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "0.02425874574963105 second"
      ],
      "text/latex": [
       "$0.02425874574963105\\ \\mathrm{second}$"
      ],
      "text/plain": [
       "0.02425874574963105 <Unit('second')>"
      ]
     },
     "execution_count": 39,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "t_final2 - t_final"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "That's not huge, which suggests we might not be able to save much time by decreasing the frontal area, or coefficient of drag, of the car."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Rolling resistance"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next we'll consider [rolling resistance](https://en.wikipedia.org/wiki/Rolling_resistance), which the force that resists the motion of the car as it rolls on tires.  The cofficient of rolling resistance, `C_rr`, is the ratio of rolling resistance to the normal force between the car and the ground (in that way it is similar to a coefficient of friction).\n",
    "\n",
    "The following function computes rolling resistance."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [],
   "source": [
    "system.set(unit_rr = 1 * N / kg)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [],
   "source": [
    "def rolling_resistance(system):\n",
    "    \"\"\"Computes force due to rolling resistance.\n",
    "    \n",
    "    system: System object\n",
    "    \n",
    "    returns: force\n",
    "    \"\"\"\n",
    "    return -system.C_rr * system.mass * system.unit_rr"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The acceleration due to rolling resistance is 0.2 (it is not a coincidence that it equals `C_rr`)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "-60.0 newton"
      ],
      "text/latex": [
       "$-60.0\\ \\mathrm{newton}$"
      ],
      "text/plain": [
       "-60.0 <Unit('newton')>"
      ]
     },
     "execution_count": 42,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "rolling_resistance(system)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "-0.2 newton/kilogram"
      ],
      "text/latex": [
       "$-0.2\\ \\frac{\\mathrm{newton}}{\\mathrm{kilogram}}$"
      ],
      "text/plain": [
       "-0.2 <Unit('newton / kilogram')>"
      ]
     },
     "execution_count": 43,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "rolling_resistance(system) / system.mass"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's a modified slope function that includes drag and rolling resistance."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [],
   "source": [
    "def slope_func3(state, t, system):\n",
    "    \"\"\"Computes the derivatives of the state variables.\n",
    "    \n",
    "    state: State object\n",
    "    t: time\n",
    "    system: System object \n",
    "    \n",
    "    returns: sequence of derivatives\n",
    "    \"\"\"\n",
    "    x, v = state\n",
    "    r_wheel, speed_ratio = system.r_wheel, system.speed_ratio\n",
    "    mass = system.mass\n",
    "    \n",
    "    omega2 = v / r_wheel * radian\n",
    "    omega1 = omega2 / speed_ratio\n",
    "    tau1 = compute_torque(omega1, system)\n",
    "    tau2 = tau1 / speed_ratio\n",
    "    F = tau2 / r_wheel\n",
    "    a_motor = F / mass\n",
    "    a_drag = drag_force(v, system) / mass\n",
    "    a_roll = rolling_resistance(system) / mass\n",
    "    \n",
    "    a = a_motor + a_drag + a_roll\n",
    "    return v, a "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And here's the run."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>values</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>success</th>\n",
       "      <td>True</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>message</th>\n",
       "      <td>A termination event occurred.</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "success                             True\n",
       "message    A termination event occurred.\n",
       "dtype: object"
      ]
     },
     "execution_count": 45,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "results3, details = run_ode_solver(system, slope_func3, events=event_func)\n",
    "details"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The final time is a little higher, but the total cost of rolling resistance is only 3/100 seconds."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "2.0646544476416953 second"
      ],
      "text/latex": [
       "$2.0646544476416953\\ \\mathrm{second}$"
      ],
      "text/plain": [
       "2.0646544476416953 <Unit('second')>"
      ]
     },
     "execution_count": 46,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "t_final3 = get_last_label(results3) * s"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "0.030300127224499374 second"
      ],
      "text/latex": [
       "$0.030300127224499374\\ \\mathrm{second}$"
      ],
      "text/plain": [
       "0.030300127224499374 <Unit('second')>"
      ]
     },
     "execution_count": 47,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "t_final3 - t_final2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So, again, there is probably not much to be gained by decreasing rolling resistance.\n",
    "\n",
    "In fact, it is hard to decrease rolling resistance without also decreasing traction, so that might not help at all."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Optimal gear ratio"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The gear ratio 13:60 is intended to maximize the acceleration of the car without causing the tires to slip.  In this section, we'll consider other gear ratios and estimate their effects on acceleration and time to reach 100 kph.\n",
    "\n",
    "Here's a function that takes a speed ratio as a parameter and returns time to reach 100 kph."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {},
   "outputs": [],
   "source": [
    "def time_to_speed(speed_ratio, params):\n",
    "    \"\"\"Computes times to reach 100 kph.\n",
    "    \n",
    "    speed_ratio: ratio of wheel speed to motor speed\n",
    "    params: Params object\n",
    "    \n",
    "    returns: time to reach 100 kph, in seconds\n",
    "    \"\"\"\n",
    "    params = Params(params, speed_ratio=speed_ratio)\n",
    "    system = make_system(params)\n",
    "    system.set(unit_rr = 1 * N / kg)\n",
    "    \n",
    "    results, details = run_ode_solver(system, slope_func3, events=event_func)\n",
    "    t_final = get_last_label(results)\n",
    "    a_initial = slope_func(system.init, 0, system)\n",
    "    return t_final"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can test it with the default ratio:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2.0646544476416953"
      ]
     },
     "execution_count": 49,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "time_to_speed(13/60, params)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can try it with different numbers of teeth on the motor gear (assuming that the axle gear has 60 teeth):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "8 1.3230554808694261\n",
      "9 1.4683740716590767\n",
      "10 1.6154033363003908\n",
      "11 1.763893473709603\n",
      "12 1.913673186217739\n",
      "13 2.0646544476416953\n",
      "14 2.216761311453768\n",
      "15 2.369962929121199\n",
      "16 2.5242340753735495\n",
      "17 2.6795453467447845\n"
     ]
    }
   ],
   "source": [
    "for teeth in linrange(8, 18):\n",
    "    print(teeth, time_to_speed(teeth/60, params))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Wow!  The speed ratio has a big effect on the results.  At first glance, it looks like we could break the world record (1.513 seconds) just by decreasing the number of teeth.\n",
    "\n",
    "But before we try it, let's see what effect that has on peak acceleration."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {},
   "outputs": [],
   "source": [
    "def initial_acceleration(speed_ratio, params):\n",
    "    \"\"\"Maximum acceleration as a function of speed ratio.\n",
    "    \n",
    "    speed_ratio: ratio of wheel speed to motor speed\n",
    "    params: Params object\n",
    "    \n",
    "    returns: peak acceleration, in m/s^2\n",
    "    \"\"\"\n",
    "    params = Params(params, speed_ratio=speed_ratio)\n",
    "    system = make_system(params)\n",
    "    a_initial = slope_func(system.init, 0, system)[1] * m/s**2\n",
    "    return a_initial"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here are the results:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "8 23.076923076923077 meter * newton / kilogram / second ** 2\n",
      "9 20.51282051282051 meter * newton / kilogram / second ** 2\n",
      "10 18.46153846153846 meter * newton / kilogram / second ** 2\n",
      "11 16.783216783216787 meter * newton / kilogram / second ** 2\n",
      "12 15.384615384615385 meter * newton / kilogram / second ** 2\n",
      "13 14.201183431952662 meter * newton / kilogram / second ** 2\n",
      "14 13.186813186813184 meter * newton / kilogram / second ** 2\n",
      "15 12.307692307692308 meter * newton / kilogram / second ** 2\n",
      "16 11.538461538461538 meter * newton / kilogram / second ** 2\n",
      "17 10.85972850678733 meter * newton / kilogram / second ** 2\n"
     ]
    }
   ],
   "source": [
    "for teeth in linrange(8, 18):\n",
    "    print(teeth, initial_acceleration(teeth/60, params))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As we decrease the speed ratio, the peak acceleration increases.  With 8 teeth on the motor gear, we could break the world record, but only if we can accelerate at 2.3 times the acceleration of gravity, which is impossible without very sticky ties and a vehicle that generates a lot of downforce."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2.354081632653061"
      ]
     },
     "execution_count": 53,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "23.07 / 9.8"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "These results suggest that the most promising way to improve the performance of the car (for this event) would be to improve traction."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
